8

我试图复制柯蒂斯在轨道力学中的一个情节,但我就是不太明白。但是,我通过切换到np.arctan2from取得了进展np.arctan

也许我实施arctan2不正确?

import pylab
import numpy as np

e = np.arange(0.0, 1.0, 0.15).reshape(-1, 1)

nu = np.linspace(0.001, 2 * np.pi - 0.001, 50000)
M2evals = (2 * np.arctan2(1, 1 / (((1 - e) / (1 + e)) ** 0.5 * np.tan(nu / 2) -
           e * (1 - e ** 2) ** 0.5 * np.sin(nu) / (1 + e * np.cos(nu)))))

fig2 = pylab.figure()
ax2 = fig2.add_subplot(111)

for Me2, _e in zip(M2evals, e.ravel()):
    ax2.plot(nu.ravel(), Me2, label = str(_e))

pylab.legend()
pylab.xlim((0, 7.75))
pylab.ylim((0, 2 * np.pi))
pylab.show()

在下图中,出现了不连续性。该函数应该是平滑的,并且在 (0, 2pi) 的 y 范围内的 0 和 2 pi 处连接,不接触 0 和 2pi。

在此处输入图像描述

教科书情节和方程式:

在此处输入图像描述

在此处输入图像描述

应 Saullo Castro 的要求,我被告知:

问题可能在于给出“原则值”作为输出的 arctan 函数。

因此,如果 x 是第二或第三象限中的角度,则 arctan(tan(x)) 不会产生 x。如果你从 x = 0 到 x = Pi 绘制 arctan(tan(x)),你会发现它在 x = Pi/2 处有一个不连续的跳跃。

对于你的情况,我相信你会写 arctan2(1, 1/arg) 而不是写 arctan(arg),其中 arg 是你的 arctan 函数的参数。这样,当 arg 变为负时,arctan2 将在第二象限而不是第四象限中产生一个角度。”

4

1 回答 1

12

通常的做法是在 的负数结果中求和 2 pi arctan()这可以有效地完成。OP 建议用 arctan2(1, 1/x) 替换 arctan(x),也由 Yay295 指出的 Maple 15 的文档提出,产生相同的结果,而无需求和 2 pi。两者如下图所示:

import pylab
import numpy as np
e = np.arange(0.0, 1.0, 0.15).reshape(-1, 1)
nu = np.linspace(0, 2*np.pi, 50000)
x =  ((1-e)/(1+e))**0.5 * np.tan(nu/2.)
x2 = e*(1-e**2)**0.5 * np.sin(nu)/(1 + e*np.cos(nu))
using_arctan = True
using_OP_arctan2 = False

if using_arctan:
    M2evals = 2*np.arctan(x) - x2
    M2evals[M2evals<0] += 2*np.pi
elif using_OP_arctan2:
    M2evals = 2 * np.arctan2(1,1/x) - x2

fig2 = pylab.figure()
ax2 = fig2.add_subplot(111)
for M2e, _e in zip(M2evals, e.ravel()):
    ax2.plot(nu.ravel(), M2e, label = str(_e))
pylab.legend(loc='upper left')
pylab.show()

在此处输入图像描述

于 2013-05-17T17:34:39.070 回答