MATLAB动态星星运转模型

MATLAB动态星星运转模型

示意图如下：

clc
close
clear
k=600;
dk=10; %平面网格密度
a=400;%%长半轴
b=200;%%半焦距
c=sqrt(a^2-b^2);%%短半轴
AA1=[600   0   b
     400,-300,250];%%质量和平面位置
xx=1.1;            %%%画图范围
bb=2;              %%动态圈数
xxx=0.6;           %%%z坐标上移比例
yyy=0.5;           %%%z坐标下移比例
S=pi()*a*c;
for i=1:361        %%生成行星轨迹及引力图
    R(i)=(b^2)/(a-c*sin(i*2*pi()/360));%%存储半径
    [fai]=jd(i,a,b);                   %%θ为椭圆上点P与焦点的连线与y轴夹角
    FF(i)=fai;                         %%存储角度
    AA1(2,2)=R(i)*cos(fai);
    AA1(2,3)=R(i)*sin(fai)-c;
    X(i)=AA1(2,2);                     %%存储空间点位X坐标
    Y(i)=AA1(2,3);                     %%存储空间点位Y坐标
    [XX,YY,Z]=SHDT(AA1 ,k ,dk);
    Data{i} = Z;
end
i=1;
[x,y,z]=sphere;     %%画星球
C=AA1(:,1)/30;      %%星球直径
figure 
set(gcf,'Unit','normalized','Position',[0,0,1,0.95])
h1=mesh(C(1)*x+AA1(i,2),C(1)*y+AA1(i,3),C(1)*z,'EdgeColor','r');%%画恒星
hold on 
h0=mesh(XX,YY,Z) ;                                              %%画引力图
h2=mesh(C(2)*x+AA1(i+1,2),C(2)*y+AA1(i+1,3),C(2)*z,'EdgeColor','g');%%画行星
% h3=plot(X,Y);
hold on
axis([-1*xx*k,xx*k,-1*xx*k,xx*k,-1*(xx-xxx)*k,(xx-yyy)*k])       %%视图范围
grid on 
title('引力模型运转系统');
view(30,30)%%%%视角
k=0;
plot(X,Y);
while k<bb             %%转的圈数
    for i=1:361
    set(h0,'XData',XX,'YData',YY,'ZData',Data{i})
    set(h1,'XData',C(1)*x+AA1(1,2),'YData',C(1)*y+AA1(1,3),'ZData',C(1)*z);
    set(h2,'XData',C(2)*x+X(i),'YData',C(2)*y+Y(i),'ZData',C(2)*z);
%     set(h3,'XData',X(i),'YData',Y(i),'ZData',[0]);
    drawnow  %现在就画
     pause(0.01)%%画图延时时间
    end
k=k+1;
end
%  close
%  plot(FF,R,'*')%%%绘制角度半径图


% clc
% close
% clear

function [XX,YY,Z]=SHDT(A1,a,da)


x=-a:da:a;
[XX,YY]=meshgrid(x);
xx=1.1;%%%画图范围
xxx=0.7;%%%z坐标下移比例
% A1=[50 0 0
%      40,-30,25];%%质量和平面位置
Z=zeros();
C=A1(:,1)/30;%%星球直径
for i=1:length(A1(:,1))
    
    AA=A1(i,:);
    [Z1]=SZ(XX,YY,AA,C(i));
    Z=Z+Z1;
end
% [x,y,z]=sphere;%%画星球
% for j=1:length(A1(:,1))
%     c=A1(j,1)/15;%%星体直径
%     mesh(c*x+A1(j,2),c*y+A1(j,3),c*z,'EdgeColor','g')
% %     hold on
% end
% Z=Z/10;
% Z=Z-max(max(Z))*ones();
% meshc(XX,YY,Z);%%画万有引力面
% axis([-1*xx*a,xx*a,-1*xx*a,xx*a,(-1*xx-xxx)*a,(xx-xxx)*a])
% view(30,30)
end

``