卫星变轨轨迹和推力模拟（单一引力源）MATLAB

 
代码阐明：

• 常量界说：界说了万有引力常数、地球和月球的质量、半径以及地月平均距离。
  
• 初始状态设置：设置卫星的初始位置、速率和姿态，以及月球的初始位置。
  
• 模拟循环：在循环中计算地球和月球对卫星的引力，模拟机动变轨推力，更新卫星的速率、位置和姿态。
  
• 绘制 3D 场景：卫星轨迹、卫星姿态和推力矢量。
  
• 辅助函数：rotationMatrix 函数用于计算旋转矩阵，drawSatellite 函数用于绘制卫星（长方体）。
  

你可以根据必要调整模拟参数和卫星的初始状态。

1. % 常量定义
2. G = 6.67430e-20; % 万有引力常数 (km^3 kg^-1 s^-2)
3. M_earth = 5.972e12; % 地球质量 (kg)
5. % 卫星初始状态
6. initial_position = [5; 0; 0]; % 初始位置 (km)，椭圆轨道长半轴的端点
7. initial_velocity = [0; 2; 0]; % 初始速度 (km/s)，垂直于长半轴方向
8. initial_attitude = [1 0 0; 0 1 0; 0 0 1]; % 初始姿态矩阵
10. % 引力源位置（假设在原点）
11. gravitational_source_position = [0; 0; 0];
13. % 模拟参数
14. total_time = 200; % 总模拟时间 (s)，根据椭圆轨道周期调整
15. dt = 0.1; % 时间步长 (s)
16. num_steps = total_time / dt;
18. % 初始化数组
19. position = zeros(3, num_steps);
20. velocity = zeros(3, num_steps);
21. attitude = zeros(3, 3, num_steps);
23. position(:, 1) = initial_position;
24. velocity(:, 1) = initial_velocity;
25. attitude(:, :, 1) = initial_attitude;
27. % 椭圆轨道参数
28. semi_major_axis = norm(initial_position); % 长半轴
29. angular_velocity = norm(initial_velocity) / semi_major_axis; % 角速度
31. % 模拟循环
32. for i = 1:num_steps - 1
33.     % 计算引力
34.     r = position(:, i) - gravitational_source_position;
35.     gravitational_force = -G * M_earth * r / norm(r)^3;
36.    
37.     % 计算当前位置的角度
38.     current_angle = atan2(position(2, i), position(1, i));
39.    
40.     % 计算向心力方向的单位向量
41.     centripetal_direction = -position(:, i) / norm(position(:, i));
42.    
43.     % 计算向心力大小（根据圆周运动公式 F = m * v^2 / r，这里假设 m = 1）
44.     centripetal_magnitude = norm(velocity(:, i))^2 / norm(position(:, i));
45.    
46.     % 计算推力
47.     thrust = centripetal_direction * centripetal_magnitude;
48.    
49.     % 计算总加速度
50.     acceleration = gravitational_force + thrust;
51.    
52.     % 更新速度和位置
53.     velocity(:, i + 1) = velocity(:, i) + acceleration * dt;
54.     position(:, i + 1) = position(:, i) + velocity(:, i) * dt;
55.    
56.     % 更新姿态（这里简单假设姿态绕推力方向旋转一个小角度）
57.     rotation_angle = 0.001; % 旋转角度 (rad)
58.     rotation_axis = thrust / norm(thrust);
59.     rotation_matrix = rotationMatrix(rotation_axis, rotation_angle);
60.     attitude(:, :, i + 1) = rotation_matrix * attitude(:, :, i);
61. end
63. % 绘制 3D 轨迹
64. figure;
65. plot3(position(1, :), position(2, :), position(3, :), 'b', 'LineWidth', 2);
66. xlabel('X (km)');
67. ylabel('Y (km)');
68. zlabel('Z (km)');
69. title('卫星变轨 3D 轨迹（考虑引力源）');
70. grid on;
72. % 绘制卫星姿态和推力矢量
73. for i = 1:10:num_steps
74.     % 绘制卫星（长方体）
75.     drawSatellite(position(:, i), attitude(:, :, i));
76.    
77.     % 确定立方体底面中心位置
78.     length = 1;
79.     width = 0.5;
80.     height = 0.2;
81.     base_center_local = [length/2; width/2; 0];
82.     base_center_global = attitude(:, :, i) * base_center_local + position(:, i);
83.    
84.     % 计算当前位置的角度
85.     current_angle = atan2(position(2, i), position(1, i));
86.    
87.     % 计算向心力方向的单位向量
88.     centripetal_direction = -position(:, i) / norm(position(:, i));
89.    
90.     % 计算向心力大小（根据圆周运动公式 F = m * v^2 / r，这里假设 m = 1）
91.     centripetal_magnitude = norm(velocity(:, i))^2 / norm(position(:, i));
92.    
93.     % 计算推力
94.     thrust = centripetal_direction * centripetal_magnitude;
95.    
96.     % 绘制推力矢量
97.     quiver3(base_center_global(1), base_center_global(2), base_center_global(3), thrust(1), thrust(2), thrust(3), 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5);
98.     % 标注矢量大小
99.     text(base_center_global(1)+thrust(1), base_center_global(2)+thrust(2), base_center_global(3)+thrust(3), num2str(centripetal_magnitude), 'Color', 'r');
100. end
102. % 旋转矩阵函数
103. function R = rotationMatrix(axis, angle)
104.     axis = axis / norm(axis);
105.     c = cos(angle);
106.     s = sin(angle);
107.     t = 1 - c;
108.     x = axis(1);
109.     y = axis(2);
110.     z = axis(3);
111.    
112.     R = [t*x^2 + c    t*x*y - s*z  t*x*z + s*y;
113.          t*x*y + s*z  t*y^2 + c    t*y*z - s*x;
114.          t*x*z - s*y  t*y*z + s*x  t*z^2 + c];
115. end
117. % 绘制卫星（长方体）函数
118. function drawSatellite(position, attitude)
119.     % 长方体尺寸
120.     length = 1;
121.     width = 0.5;
122.     height = 0.2;
123.    
124.     % 长方体顶点坐标
125.     vertices = [0 0 0;
126.                 length 0 0;
127.                 length width 0;
128.                 0 width 0;
129.                 0 0 height;
130.                 length 0 height;
131.                 length width height;
132.                 0 width height];
133.    
134.     % 应用姿态变换
135.     vertices = vertices * attitude';
136.    
137.     % 平移到当前位置
138.     vertices = bsxfun(@plus, vertices, position');
139.    
140.     % 绘制长方体
141.     faces = [1 2 3 4; 5 6 7 8; 1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8];
142.     patch('Vertices', vertices, 'Faces', faces, 'FaceColor', 'g', 'EdgeColor', 'k');
143.     hold on;
144. end   

复制代码

 

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。