韦伯望远镜的拉格朗日点计算推导过程，包含MATLAB和python运动轨迹仿真代码

研究过程

起源与提出：1687 年牛顿提出 “三体问题”，旨在研究三个可视为质点的天体在相互之间万有引力作用下的运动规律，但因运动方程过于复杂，难以得到完全解。
欧拉的贡献1：1767 年，瑞士数学家莱昂哈德・欧拉根据旋转的二体引力场推算出了其中的三个特解，即 L1、L2 和 L3。
拉格朗日的突破1：1772 年，法国数学家约瑟夫・路易・拉格朗日推算出了剩下的两个特解 L4 和 L5。至此，平面圆型限制性三体问题的五个特解全部被找到，后来人们将这五个点统称为拉格朗日点。
观测证实1：1906 年，天文学家马克斯・沃尔夫发现一颗位于火星和木星间主带以外的小行星，它与木星、太阳构成等边三角形，位于木星前方的 L4 点。同年还发现 617 号小行星在木星后方的 L5 点。20 世纪 80 年代，在土星和它的大卫星构成的运动系统中也发现类似等边三角形。截至 2009 年，天文学家在木星的 L4 和 L5 周围各发现超过 1000 颗小行星，进一步证实了拉格朗日点的存在。

重要意义

在天文学领域
- 解释天体分布：有助于解释为何某些小行星和卫星能稳定存在于行星或太阳的特定轨道上，如木星的特洛伊小行星群位于木星的 L4 和 L5 点。
- 观测优势1：日地系统的 L1 点能提供不间断的太阳视野，是观测太阳的绝佳位置；L2 点有稳定的热力学环境以及与地球相对固定的构型，利于开展天文观测任务，如詹姆斯・韦伯空间望远镜就位于日地系统的 L2 点。
在航天领域
- 节省能源：航天器位于拉格朗日点时，受两个大天体引力作用，可保持相对稳定的轨道，能以最小的能量消耗维持轨道，节省大量燃料，降低航天任务成本。
- 作为星际探测中转站1：可作为星际探测的理想中转站，如月球 - 地球系统的 L1 点可作为地月转移的中继点，为航天器的轨道转移和姿态调整等提供便利，有助于更深入地开展深空探测任务。
在天体力学理论方面
- 推动三体问题研究：是限制性三体问题的重要特解，为三体问题的研究提供了特殊情况和突破口，帮助科学家更好地理解和处理复杂的天体力学问题，推动了天体力学理论的发展。
- 验证引力理论：对拉格朗日点的研究和观测，是对牛顿万有引力定律等引力理论的有效验证，进一步巩固了引力理论在天体力学中的地位，也为广义相对论等更深入的引力理论研究提供了实践基础。
在未来太空开发方面
- 潜在太空殖民地选址：拉格朗日点尤其是 L5 点相对稳定，可作为未来太空殖民地的潜在位置，为人类在太空中长期居住和建设提供可能，为解决地球资源、人口等问题提供新的思路和方向。
- 资源开发的战略点：未来若对小行星等天体资源进行开发，拉格朗日点可作为资源运输、存储和加工的战略节点，便于对资源进行集中处理和利用，为太空资源开发产业的发展提供支撑。

1、推导过程

3. Python 代码实现拉格朗日点计算和运动轨迹仿真

python">import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import newton# 定义常量
G = 6.67430e-11  # 万有引力常数 (m^3 kg^-1 s^-2)
M_sun = 1.989e30  # 太阳质量 (kg)
M_earth = 5.972e24  # 地球质量 (kg)
a = 1.496e11  # 地球到太阳的平均距离 (m)
T = 365.25 * 24 * 3600  # 地球绕太阳的公转周期 (s)
omega = 2 * np.pi / T  # 地球绕太阳的公转角速度 (rad/s)
mu = M_earth / (M_sun + M_earth)# 定义求解 L2 点的方程
def f(x):return -omega**2 * x + G * M_sun / (x + mu * a)**2 + G * M_earth / (x - (1 - mu) * a)**2# 使用牛顿 - 拉夫逊方法求解 L2 点
L2 = newton(f, a + 1e9)print(f"L2 点到太阳的距离: {L2:.2e} m")# 运动轨迹仿真
t = np.linspace(0, T, 1000)
x_sun = -mu * a * np.cos(omega * t)
y_sun = -mu * a * np.sin(omega * t)
x_earth = (1 - mu) * a * np.cos(omega * t)
y_earth = (1 - mu) * a * np.sin(omega * t)
x_L2 = L2 * np.cos(omega * t)
y_L2 = L2 * np.sin(omega * t)# 绘制运动轨迹
plt.figure(figsize=(10, 10))
plt.plot(x_sun, y_sun, label='Sun')
plt.plot(x_earth, y_earth, label='Earth')
plt.plot(x_L2, y_L2, label='L2 Point')
plt.scatter(0, 0, color='black', label='Barycenter')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.title('Sun - Earth - L2 Point Orbits')
plt.legend()
plt.axis('equal')
plt.grid(True)
plt.show()

效果如下

代码解释

常量定义：定义了万有引力常数、太阳质量、地球质量、地球到太阳的平均距离、地球绕太阳的公转周期和公转角速度等常量。
求解 L2 点：定义了求解 L2 点的方程 f(x)，并使用 scipy.optimize.newton 函数求解该方程。
运动轨迹仿真：计算太阳、地球和 L2 点在一段时间内的位置，并使用 matplotlib 库绘制它们的运动轨迹。

注意事项

上述代码假设韦伯望远镜精确位于 L2 点，实际情况中，望远镜会围绕 L2 点进行小幅度的轨道保持操作。
代码中的数值计算是基于简化的二体问题模型，实际的太阳系是一个多体系统，需要考虑更多的因素。

4. MATLAB 代码实现拉格朗日点计算和运动轨迹仿真

代码如下

G = 6.67430e-11; % 万有引力常数 (m^3 kg^-1 s^-2)
M_sun = 1.989e30; % 太阳质量 (kg)
M_earth = 5.972e24; % 地球质量 (kg)
a = 1.496e11; % 地球到太阳的平均距离 (m)
T = 365.25 * 24 * 3600; % 地球绕太阳的公转周期 (s)
omega = 2 * pi / T; % 地球绕太阳的公转角速度 (rad/s)
mu = M_earth / (M_sun + M_earth);% 定义求解 L2 点的函数
f = @(x) -omega^2 * x + G * M_sun / (x + mu * a)^2 + G * M_earth / (x - (1 - mu) * a)^2;% 使用牛顿 - 拉夫逊方法求解 L2 点
x0 = a + 1e9; % 初始猜测值
L2 = fzero(f, x0); 
% 创建 3D 图形窗口
figure('Position', [100, 100, 800, 600]);
ax = gca;
ax.XLim = [-3, 3];
ax.YLim = [-3, 3];
ax.ZLim = [-1, 1];
xlabel('X');
ylabel('Y');
zlabel('Z');
title('3D Celestial Bodies Animation');
grid on;% 初始化点
x = 0; y = 0; z = 0;
h = plot3(x, y, z, 'ro'); % 'ro' 表示红色圆点
hold on
x2 = 0; y2 = 0; z2 = 0;
m = plot3(x2, y2, z2, 'bo'); 
hold on
x3 = 0; y3 = 0; z3 = 0;
n = plot3(x3, y3, z3, 'yo'); 
hold on
legend('h', 'm', 'n');
grid on% 动画循环% 计算不同时间点的位置
t = linspace(0, T, 1000);
x_sun = -mu * a * cos(omega * t);
y_sun = -mu * a * sin(omega * t);
z_sun = zeros(size(t));x_earth = (1 - mu) * a * cos(omega * t);
y_earth = (1 - mu) * a * sin(omega * t);
z_earth = zeros(size(t));x_L2 = L2 * cos(omega * t);
y_L2 = L2 * sin(omega * t);
z_L2 = zeros(size(t));set(h, 'XData', x_sun, 'YData', y_sun, 'ZData', z_sun); % 更新点的位置
set(m, 'XData', x_earth, 'YData', y_earth, 'ZData', z_earth); % 更新点的位置
set(n, 'XData', x_L2, 'YData', y_L2, 'ZData', z_L2); % 更新点的位置drawnow; % 更新图形
pause(0.1); % 暂停一段时间以观察动画效果

效果如下

代码解释：

常量定义：定义了仿真所需的各种物理常量，如万有引力常数、太阳和地球的质量、日地距离、公转周期等。
求解 L2 点：定义了一个匿名函数 f 来表示求解 L2 点的方程，然后使用 fzero 函数求解该方程得到 L2 点的位置。
位置计算：计算太阳、地球和 L2 点在不同时间点的三维位置，存储在对应的 x、y、z 数组中。这里假设所有天体都在同一平面（z 坐标为 0）上运动。
3D 图形初始化：创建一个 3D 图形窗口，设置坐标轴范围、标签和标题，并初始化表示太阳、地球和 L2 点的图形元素。
动画循环：使用 for 循环遍历所有时间点，在每一帧中更新物体的位置，并使用 drawnow limitrate 函数刷新图形，从而实现动画效果。