基于MATLAB实现的LISA轨道运动仿真+使用说明文档.rar

function [Phi, E, Sp1path, Sp2path, Sp3path, C_total, Spsup, Spinf, R] = simu_LISA_orbits() R = 15; % *10^6 km (UA) 到太阳距离Sun E = [R;0;0]; % 初始地球位置= [x_T(1);y_T(1);z_T(1)] % 3架LISA航天器的初始位置 Delta_phi = -pi/6; L = 5; % 航天器之间的距离 % 旋转矩阵 (E,x,y,z) to (cdm,x,y,z) Mr0 = [cos(Delta_phi),-sin(Delta_phi),0; sin(Delta_phi),cos(Delta_phi),0; 0,0,1]; t = (1:360)'; % 地球轨道一周 Omega = 2*pi/t(end); Phi = 0:Omega:2*pi; % 圆形路径/轨迹建模 % 地球 x_E = R.*cos(Phi); % 逆时针旋转方式 y_E = R.*sin(Phi); z_E = zeros(1,length(x_E)); E = [x_E; y_E; z_E]; %(cdm,x-y-z) plan to (S,x-y-z) plan t_v = [R;0;0]; % LISA的旋转矩阵(cdm,1,2,3) to (cdm,x-y-z) plan Mr2 = [cos(pi/3),0,-sin(pi/3); 0,1,0; sin(pi/3),0,cos(pi/3)]; % 从一个LISA航天器到其后续连续的旋转矩阵 Mr3 = [cos(-2*pi/3),-sin(-2*pi/3),0; sin(-2*pi/3),cos(-2*pi/3),0; 0,0,1]; % 初始航天器位置 Sp1 = [L/sqrt(3);0; 0]; % 航天器轨迹 for k = 1:length(E) Phi_k = Phi(1,k); % 属于LISA圆的点的旋转矩阵 Mr5 = [cos(Phi_k),-sin(Phi_k),0 sin(Phi_k),cos(Phi_k),0 0,0,1]; % 一周LISA三角形的旋转矩阵 Mr6 = Mr5.*[1,-1,1;-1,1,1;1,1,1]; Sp1path(:,k) = Mr0*(E(:,k)+Mr5*Mr2*Mr6*Sp1); end Sp2path = Mr3\Sp1path; Sp3path = Mr3*Sp1path; % 开始 Sp2path = [Sp2path(:,2*length(Sp2path)/3+1:end) Sp2path(:,1:2*length(Sp2path)/3)]; Sp3path = [Sp3path(:,length(Sp3path)/3+1:end) Sp3path(:,1:length(Sp3path)/3)]; % 连接点（调整大小正常） Phi = [Phi 2*pi]; Sp1path = [Sp1path Sp1path(:,1)]; Sp2path = [Sp2path Sp2path(:,1)]; Sp3path = [Sp3path Sp3path(:,1)]; E = [E E(:,1)]; x_E = E(1,:); y_E = E(2,:); z_E = E(3,:); % LISA 轨迹圈 C = zeros(3,length(Phi)); for k = 1:length(Phi) Phi_k = Phi(1,k); % LISA圆的点的旋转矩阵 Mr5 = [cos(Phi_k),-sin(Phi_k),0 sin(Phi_k),cos(Phi_k),0 0,0,1]; % 一周期LISA三角形的旋转矩阵 Mr6 = Mr5.*[1,-1,1;-1,1,1;1,1,1]; C(:,k) = Mr0*(t_v+Mr2*Mr6*Sp1); end C_total = zeros(3,length(Phi),length(Phi)); C_total(:,:,1) = C(:,:); % 初始圆 for n = 2:length(Phi) Phi_n = Phi(1,n); % 属于LISA圆的点的旋转矩阵 Mr5 = [cos(Phi_n),-sin(Phi_n),0 sin(Phi_n),cos(Phi_n),0 0,0,1]; C_total(:,:,n) = Mr5*C_total(:,:,1); end % LISA轨道路径限制 x_Spsup = x_E+(L*cos(pi/3)/sqrt(3))*cos(Phi); y_Spsup = y_E+(L*cos(pi/3)/sqrt(3))*sin(Phi); z_Spsup = z_E+L*sin(pi/3)/sqrt(3); Spsup = [x_Spsup;y_Spsup;z_Spsup]; Spsup = Mr0*Spsup; x_Spinf = x_E-(L*cos(pi/3)/sqrt(3))*cos(Phi); y_Spinf = y_E-(L*cos(pi/3)/sqrt(3))*sin(Phi); z_Spinf = z_E-L*sin(pi/3)/sqrt(3); Spinf = [x_Spinf;y_Spinf;z_Spinf]; Spinf = Mr0*Spinf; end