% 水星距离太阳最远处的距离 a = 0.6982e11; % 单位：m % 线速度 v = 3.886e4; % 单位：m/s % 椭圆的短轴长度 b = a * v; c=sqrt(a^2-b^2); % Kepler 第三定律计算周期 G = 6.67430e-11; % 引力常数，单位：m^3/(kg*s^2) M = 1.989e30; % 太阳的质量，单位：kg T = 2 * pi * sqrt(a^3 / (G * M)); % 单位：s % 水星到太阳的最近距离 closest_distance = b; % 第 50 天结束时水星的位置 time = 50 * 24 * 60 * 60; % 50天的秒数 mean_angular_velocity = 2 * pi / T; % 平均角速度 angle = mean_angular_velocity * time; % 角度 x = a * cos(angle); y = b * sin(angle); % 绘制水星绕太阳运行的轨道曲线 t = linspace(0, 2 * pi, 1000); orbit_x = a * cos(t); orbit_y = b * sin(t); % 绘图 figure; plot(orbit_x, orbit_y, 'b'); hold on; plot(x, -y, 'ro'); xlabel('x'); ylabel('y'); title('Orbit of Mercury'); legend('Orbit', 'Mercury Position');，在此代码上改进一点，使最后的图形会打出椭圆轨道的两个焦点

要画出椭圆轨道的两个焦点，需要先求出椭圆的长轴和焦距。长轴可以直接用距离太阳最远处的距离 a 计算得到。而焦距可以用勾股定理算出：

f = sqrt(a^2 - b^2);

接下来，在绘图前加上以下代码，用红色的点标出两个焦点：

focus1 = [-f, 0];
focus2 = [f, 0];
plot(focus1(1), focus1(2), 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r');
plot(focus2(1), focus2(2), 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r');

最终的代码如下：

% 水星距离太阳最远处的距离
a = 0.6982e11; % 单位：m

% 线速度
v = 3.886e4; % 单位：m/s

% 椭圆的短轴长度
b = a * v; 

% 椭圆的焦距
f = sqrt(a^2 - b^2);

% Kepler 第三定律计算周期
G = 6.67430e-11; % 引力常数，单位：m^3/(kg*s^2)
M = 1.989e30; % 太阳的质量，单位：kg
T = 2 * pi * sqrt(a^3 / (G * M)); % 单位：s

% 水星到太阳的最近距离
closest_distance = b;

% 第 50 天结束时水星的位置
time = 50 * 24 * 60 * 60; % 50天的秒数
mean_angular_velocity = 2 * pi / T; % 平均角速度
angle = mean_angular_velocity * time; % 角度
x = a * cos(angle);
y = b * sin(angle);

% 绘制水星绕太阳运行的轨道曲线
t = linspace(0, 2 * pi, 1000);
orbit_x = a * cos(t);
orbit_y = b * sin(t);

% 绘图
figure;
plot(orbit_x, orbit_y, 'b');
hold on;
plot(x, -y, 'ro');

% 绘制焦点
focus1 = [-f, 0];
focus2 = [f, 0];
plot(focus1(1), focus1(2), 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r');
plot(focus2(1), focus2(2), 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r');

xlabel('x');
ylabel('y');
title('Orbit of Mercury');
legend('Orbit', 'Mercury Position', 'Foci');