% 初始参数  
m_satellite = 500; % 卫星质量 (kg)  
m_rocket = 100;    % 火箭壳体质量 (kg)  
v_initial = 7000;  % 初始共同速度 (m/s)  
v_rel = 1800;      % 相对速度 (m/s)  
  
% 根据动量守恒计算分离后的速度  
v_satellite = (m_satellite*v_initial + m_rocket*(v_initial + v_rel)) / (m_satellite + m_rocket);  
v_rocket = v_satellite - v_rel;  
  
% 输出结果  
fprintf('分离后卫星的速度: %.2f m/s\n', v_satellite);  
fprintf('分离后火箭壳体的速度: %.2f m/s\n', v_rocket);  
  
% 假设它们仍在圆形轨道上运动（实际情况可能不同）  
% 可以计算新的轨道半径（仅用于示意，实际情况会更复杂）  
G = 6.67430e-11; % 万有引力常数 (m^3 kg^-1 s^-2)  
M_earth = 5.972e24; % 地球质量 (kg)  
r_initial = (G * M_earth) / (v_initial^2); % 初始轨道半径  
r_satellite = (G * M_earth) / (v_satellite^2); % 卫星新轨道半径  
r_rocket = (G * M_earth) / (v_rocket^2); % 火箭壳体新轨道半径  
  
% 输出轨道半径  
fprintf('卫星的新轨道半径: %.2e m\n', r_satellite);  
fprintf('火箭壳体的新轨道半径: %.2e m\n', r_rocket);
/*
代码中计算的“新轨道半径”是基于简化的物理假设，即物体在圆周运动中，其轨道半径 r 与速度 v 的关系是通过万有引力提供向心力这一公式得出的。具体来说，对于一个在地球附近做圆周运动的物体，其轨道半径 r 和速度 v 满足以下关系：

[ r = \frac{GM}{v^2} ]

其中 G 是万有引力常数，M 是地球的质量。这个公式表明，速度越大，轨道半径越小，这是在假设物体仍然保持在圆周轨道上的情况下得出的结论。

然而，这个结论在实际太空环境中并不总是准确的，因为它忽略了其他可能影响轨道的因素，比如大气阻力（在低地球轨道中）、其他天体的引力摄动，以及轨道变化时可能出现的椭圆轨道效应。

在你的代码中，我们根据动量守恒计算了卫星和火箭壳体分离后的新速度。由于卫星的速度增加了，根据上面的公式，计算出的“新轨道半径”就会减小。但这并不意味着卫星真的会进入一个更小的圆轨道，而是说如果卫星仍然保持在圆周轨道上运动，那么它需要的向心力将会对应一个更小的半径。

实际上，当卫星的速度发生变化时，它更可能进入一个椭圆轨道，而不是一个更小的圆轨道。在这个椭圆轨道上，近地点可能比原来的圆轨道半径要小，而远地点则可能要比原来的轨道半径大。要准确描述这种轨道变化，需要使用更复杂的轨道动力学模型。

因此，代码中的计算结果应当被视为一种理论上的、简化的模型预测，它可能并不反映真实世界中卫星轨道变化的复杂性。在实际应用中，航天工程师会使用精确的轨道分析工具和模型来预测和规划卫星的轨道。
*/

clear all

% 初始参数  
global G M_earth; 
G = 6.67430e-11; % 万有引力常数 (N*(m^2)/(kg^2))  
M_earth = 5.972e24; % 地球质量 (kg)  
m_satellite = 500; % 卫星质量 (kg)，在这个问题中其实不重要，因为不影响轨迹  
 m_rocket = 100; %火箭壳体的质量是100 kg
   
% 初始轨道参数  
v_initial = 7000; % 初始的速度是 7000 m/s 
r_initial = G*M_earth/ v_initial^2 ; % 圆周运动计算  F=mv^2/r , F=GMm/r^2 联立求解

v_rel = 1800;      % 相对速度 (m/s) 

% 根据动量守恒计算分离后的速度  
v_satellite = (m_satellite*v_initial + m_rocket*(v_initial + v_rel)) / (m_satellite + m_rocket);  
v_rocket = v_satellite - v_rel;  
  

% 设置时间范围和初始条件  
tspan = [0, 3600*3]; % 模拟3小时的时间，可以根据需要调整  
y0 = [r_initial; 0; 0; 0; v_satellite; 0]; % 初始位置在x轴上，初始的速度假设在y轴上，这样速度与半径垂直  
  
% 使用ode45求解常微分方程  
[t, y] = ode45(@satellite_dynamics, tspan, y0);  
  
% 提取位置和速度数据  
r_x = y(:,1);  
r_y = y(:,2);  
r_z = y(:,3);  
  
% 绘制卫星运动轨迹  
figure;  
plot3(r_x, r_y, r_z, 'b-');  
xlabel('X (m)');  
ylabel('Y (m)');  
zlabel('Z(m)'); 
title('Satellite Orbit after Velocity Change');  
text(r_x(1),r_y(1),r_z(1), '起点');

% 绘制卫星到地球的距离随时间变化的图  

distances = sqrt(r_x.^2 + r_y.^2 + r_z.^2);  
figure;  
plot(t, distances, 'r-');  
xlabel('Time (s)');  
ylabel('Distance from Earth (m)');  
title('Distance of Satellite from Earth over Time');  


% 定义运动方程  
function dy = satellite_dynamics(t, y)  
	global  G M_earth;
    r = y(1:3); % 卫星位置向量  
    v = y(4:6); % 卫星速度向量  
    r_magnitude = norm(r); % 到地心的距离  
    acceleration = -G * M_earth * r / r_magnitude^3; % 万有引力加速度  
    dy = [v; acceleration]; % 返回速度和加速度  
end