卫星变轨轨迹和推力模拟（单一引力源）MATLAB
代码说明：
- 常量定义：定义了万有引力常数、地球和月球的质量、半径以及地月平均距离。
- 初始状态设置：设置卫星的初始位置、速度和姿态，以及月球的初始位置。
- 模拟循环：在循环中计算地球和月球对卫星的引力，模拟机动变轨推力，更新卫星的速度、位置和姿态。
- 绘制 3D 场景：卫星轨迹、卫星姿态和推力矢量。
- 辅助函数：
rotationMatrix函数用于计算旋转矩阵，drawSatellite函数用于绘制卫星（长方体）。
你可以根据需要调整模拟参数和卫星的初始状态。
% 常量定义
G = 6.67430e-20; % 万有引力常数 (km^3 kg^-1 s^-2)
M_earth = 5.972e12; % 地球质量 (kg)% 卫星初始状态
initial_position = [5; 0; 0]; % 初始位置 (km)，椭圆轨道长半轴的端点
initial_velocity = [0; 2; 0]; % 初始速度 (km/s)，垂直于长半轴方向
initial_attitude = [1 0 0; 0 1 0; 0 0 1]; % 初始姿态矩阵% 引力源位置（假设在原点）
gravitational_source_position = [0; 0; 0];% 模拟参数
total_time = 200; % 总模拟时间 (s)，根据椭圆轨道周期调整
dt = 0.1; % 时间步长 (s)
num_steps = total_time / dt;% 初始化数组
position = zeros(3, num_steps);
velocity = zeros(3, num_steps);
attitude = zeros(3, 3, num_steps);position(:, 1) = initial_position;
velocity(:, 1) = initial_velocity;
attitude(:, :, 1) = initial_attitude;% 椭圆轨道参数
semi_major_axis = norm(initial_position); % 长半轴
angular_velocity = norm(initial_velocity) / semi_major_axis; % 角速度% 模拟循环
for i = 1:num_steps - 1% 计算引力r = position(:, i) - gravitational_source_position;gravitational_force = -G * M_earth * r / norm(r)^3;% 计算当前位置的角度current_angle = atan2(position(2, i), position(1, i));% 计算向心力方向的单位向量centripetal_direction = -position(:, i) / norm(position(:, i));% 计算向心力大小（根据圆周运动公式 F = m * v^2 / r，这里假设 m = 1）centripetal_magnitude = norm(velocity(:, i))^2 / norm(position(:, i));% 计算推力thrust = centripetal_direction * centripetal_magnitude;% 计算总加速度acceleration = gravitational_force + thrust;% 更新速度和位置velocity(:, i + 1) = velocity(:, i) + acceleration * dt;position(:, i + 1) = position(:, i) + velocity(:, i) * dt;% 更新姿态（这里简单假设姿态绕推力方向旋转一个小角度）rotation_angle = 0.001; % 旋转角度 (rad)rotation_axis = thrust / norm(thrust);rotation_matrix = rotationMatrix(rotation_axis, rotation_angle);attitude(:, :, i + 1) = rotation_matrix * attitude(:, :, i);
end% 绘制 3D 轨迹
figure;
plot3(position(1, :), position(2, :), position(3, :), 'b', 'LineWidth', 2);
xlabel('X (km)');
ylabel('Y (km)');
zlabel('Z (km)');
title('卫星变轨 3D 轨迹（考虑引力源）');
grid on;% 绘制卫星姿态和推力矢量
for i = 1:10:num_steps% 绘制卫星（长方体）drawSatellite(position(:, i), attitude(:, :, i));% 确定立方体底面中心位置length = 1;width = 0.5;height = 0.2;base_center_local = [length/2; width/2; 0];base_center_global = attitude(:, :, i) * base_center_local + position(:, i);% 计算当前位置的角度current_angle = atan2(position(2, i), position(1, i));% 计算向心力方向的单位向量centripetal_direction = -position(:, i) / norm(position(:, i));% 计算向心力大小（根据圆周运动公式 F = m * v^2 / r，这里假设 m = 1）centripetal_magnitude = norm(velocity(:, i))^2 / norm(position(:, i));% 计算推力thrust = centripetal_direction * centripetal_magnitude;% 绘制推力矢量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);% 标注矢量大小text(base_center_global(1)+thrust(1), base_center_global(2)+thrust(2), base_center_global(3)+thrust(3), num2str(centripetal_magnitude), 'Color', 'r');
end% 旋转矩阵函数
function R = rotationMatrix(axis, angle)axis = axis / norm(axis);c = cos(angle);s = sin(angle);t = 1 - c;x = axis(1);y = axis(2);z = axis(3);R = [t*x^2 + c t*x*y - s*z t*x*z + s*y;t*x*y + s*z t*y^2 + c t*y*z - s*x;t*x*z - s*y t*y*z + s*x t*z^2 + c];
end% 绘制卫星（长方体）函数
function drawSatellite(position, attitude)% 长方体尺寸length = 1;width = 0.5;height = 0.2;% 长方体顶点坐标vertices = [0 0 0;length 0 0;length width 0;0 width 0;0 0 height;length 0 height;length width height;0 width height];% 应用姿态变换vertices = vertices * attitude';% 平移到当前位置vertices = bsxfun(@plus, vertices, position');% 绘制长方体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];patch('Vertices', vertices, 'Faces', faces, 'FaceColor', 'g', 'EdgeColor', 'k');hold on;
end