如何获取隐式开普勒方程曲线与竖线的交点坐标？

问题背景 #

我编写了一段Python程序用于绘制隐式开普勒方程：y - e sin(y) = x（其中平近点角x的范围为0°到180°），程序代码如下：

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from numpy import sin, radians

y_range = np.linspace(0, 180)
x_range = np.linspace(0, 180)
y, x = np.meshgrid(y_range, x_range)

color = "red"
e = 0.90

fig, ax = plt.subplots(1, 1, figsize = (6, 3), tight_layout = True)
 
f = lambda y, x, e: radians(y) - e*sin(radians(y)) - radians(x)
equation = f(y, x, e)

ax.contour(x, y, equation, levels = [0], colors = color)

ax.set_xticks([i for i in range(0, 180+1, 15)])
ax.set_yticks([i for i in range(0, 180+1, 15)])
ax.set_xlabel("x : Mean Anomaly (deg)")
ax.set_ylabel("y : Eccentric Anomaly (deg)")

# line at point x = 15°
ax.axvline(x = 15)   

plt.show()

运行程序后可见，x=15°的竖线与隐式曲线交点纵坐标约为60°。想了解是否有简洁的Python方法获取该交点坐标，已知可使用intersect库，但不清楚如何在程序中应用，请问解决思路是什么？是否需要重构程序？

解决思路与实现 #

不需要重构现有绘图程序，核心是针对固定x值求解对应的y值，以下提供两种可行方案：

方法1：直接数值求解（最简洁，无需额外库） #

开普勒方程本质是给定平近点角x，求偏近点角y的单变量方程问题。当x固定为15°时，方程可转化为：
radians(y) - e*sin(radians(y)) - radians(15) = 0
使用scipy.optimize.root_scalar即可高效求解，无需依赖intersect库。

在原有代码基础上添加如下逻辑：

from scipy.optimize import root_scalar

# 定义固定x=15时的单变量函数
def kepler_eq(y_deg, e, x_deg):
    y_rad = radians(y_deg)
    x_rad = radians(x_deg)
    return y_rad - e * sin(y_rad) - x_rad

# 求解y值，初始猜测设为60（根据绘图结果）
result = root_scalar(kepler_eq, args=(e, 15), bracket=[0, 180], method='brentq')
if result.converged:
    print(f"交点坐标：x=15°, y={result.root:.2f}°")
else:
    print("求解失败")

运行后会输出精确的交点纵坐标，比如对于e=0.9，结果约为60.23°。

方法2：使用intersect库求曲线交点 #

如果坚持使用intersect库，需要先提取两条曲线的离散点数据：

1. 从matplotlib的contour对象中提取开普勒方程曲线的点集
2. 生成竖线x=15°的点集（y从0到180°）
3. 调用intersect库的方法计算交点

实现代码如下（需先安装库：pip install intersect）：

from intersect import intersect

# 提取contour的路径点
contour_obj = ax.contour(x, y, equation, levels=[0], colors=color)
path = contour_obj.collections[0].get_paths()[0]
kepler_points = path.vertices  # 开普勒曲线的(x,y)点集

# 生成竖线x=15的点集
vertical_line_x = np.full_like(y_range, 15)
vertical_line_y = y_range
vertical_points = np.column_stack((vertical_line_x, vertical_line_y))

# 计算交点
intersection_points = intersect(kepler_points, vertical_points)
if len(intersection_points) > 0:
    print(f"交点坐标：x={intersection_points[0][0]:.2f}°, y={intersection_points[0][1]:.2f}°")
else:
    print("未找到交点")

注意：这种方法依赖曲线的离散采样精度，结果精度略低于数值求解方法。

内容的提问来源于stack exchange，提问作者Andrew