c# - 从状态向量中找到真正的异常

Question

我正在尝试将状态向量(位置和速度)转换为开普勒元素，但是我遇到了问题，在尝试计算真正的异常时，负速度或位置会给我错误的结果。

以下是我尝试计算真实异常的不同方法:

/// <summary>
/// https://en.wikipedia.org/wiki/True_anomaly#From_state_vectors
/// </summary>
public static double TrueAnomaly(Vector4 eccentVector, Vector4 position, Vector4 velocity)
{
    var dotEccPos = Vector4.Dot(eccentVector, position);
    var talen = eccentVector.Length() * position.Length();
    talen = dotEccPos / talen;
    talen = GMath.Clamp(talen, -1, 1);
    var trueAnomoly = Math.Acos(talen);

    if (Vector4.Dot(position, velocity) < 0)
        trueAnomoly = Math.PI * 2 - trueAnomoly;

    return trueAnomoly;
}

//sgp = standard gravitational parameter
public static double TrueAnomaly(double sgp, Vector4 position, Vector4 velocity)
{
    var H = Vector4.Cross(position, velocity).Length();
    var R = position.Length();
    var q = Vector4.Dot(position, velocity);  // dot product of r*v
    var TAx = H * H / (R * sgp) - 1;
    var TAy = H * q / (R * sgp);
    var TA = Math.Atan2(TAy, TAx);
    return TA;
}

public static double TrueAnomalyFromEccentricAnomaly(double eccentricity, double eccentricAnomaly)
{
    var x = Math.Sqrt(1 - Math.Pow(eccentricity, 2)) * Math.Sin(eccentricAnomaly);
    var y = Math.Cos(eccentricAnomaly) - eccentricity;
    return Math.Atan2(x, y);
}

public static double TrueAnomalyFromEccentricAnomaly2(double eccentricity, double eccentricAnomaly)
{
    var x = Math.Cos(eccentricAnomaly) - eccentricity;
    var y = 1 - eccentricity * Math.Cos(eccentricAnomaly);
    return Math.Acos(x / y);
}

编辑:Spectre 指出的另一种方法:

public static double TrueAnomaly(Vector4 position, double loP)
{
    return Math.Atan2(position.Y, position.X) - loP; 
}

位置都是相对于父实体的。

如果 position.x、position.y 和 velocity.y 都为正，则这些函数都一致。如何解决这些问题，以便在位置和速度为负时获得一致的结果？

澄清一下:我的角度似乎有点正确，只是根据位置和/或速度矢量指向错误的象限。

<删除>是的所以我错了，毕竟以上所有都返回正确的值。

所以我发现了一个边缘情况，上面的大部分计算都失败了。给定位置和速度:

pos = new Vector4() { X = -0.208994076275941, Y = 0.955838328099748 };
vel = new Vector4() { X = -2.1678187689294E-07, Y = -7.93096769486992E-08 };

我得到一些奇怪的结果，即 ~ -31.1 度，当我认为它应该返回 `31.1(非负)。其中一个返回 ~ 328.8。
然而用这个位置和速度测试结果似乎没问题:

pos = new Vector4() { X = -0.25, Y = 0.25 };
vel = new Vector4() { X = Distance.KmToAU(-25), Y = Distance.KmToAU(-25) };

请参阅我的答案以获取有关我如何测试以及我对其他一些变量使用的数学的额外代码。

我在这个问题上兜圈子。这是我现有代码中的一个错误的结果，该错误在某些情况下出现，但在其他情况下不会出现。我想现在真正的问题是为什么我得到的位置/速度高于我的预期或彼此不匹配的不同结果？

Answer 1

假设 2D 情况...我的做法不同:

1. 计算半轴半径和旋转

   因此您需要记住整个轨道并在其上找到 2 个最远的点，即长轴 a。短轴 b 通常与长轴成 90 度角，但要确保只有 2 个垂直于您的轨道上与长轴最远的点。所以现在你有两个半轴。初始旋转由 atan2 从主轴计算。

2. 计算真实异常 E

   所以如果中心是 x0,y0(a,b 的交点或两者的中心点)初始旋转是 ang0(角度a) 并且你在轨道上的点是 x,y 那么:

   E = atan2(y-y0,x-x0) - ang0
   

然而，为了将牛顿/达朗贝尔物理学与开普勒轨道参数相匹配，您需要像我在这里所做的那样提高积分精度:

• Is it possible to make realistic n-body solar system simulation in matter of size and mass?

请参阅那里的[Edit3] Improving Newton D'ALembert integration precision even more。

有关更多信息和方程式，请参阅:

• Solving Kepler's equation

[Edit1] 所以你想计算 V 我是这样看的:

当您获得相对于父级的坐标时，您可以假设它们已经位于焦点中心，因此不再需要 x0,y0。如果你想要高精度并且有超过 2 个物体(焦点质量 + 物体 + 接近物体，如卫星)，那么母质量将不再位于轨道的焦点但靠近它......并进行补救你需要使用真实的焦点位置，所以 x0,y0 再次......那么怎么做:

1. 计算中心点(cx,cy)和a,b半轴

   所以它与前面的文本相同。

2. 在轨道轴对齐坐标中计算焦点(x0,y0)

   简单:

   x0 = cx + sqrt( a^2 + b^2 );
   y0 = cy;
   

3. a的初始角度ang0

   设xa,ya为轨道与长轴a在速度较大的一侧(靠近父对象焦点)的交点。然后:

   ang0 = atan2( ya-cy , xa-cx );
   

4. 最后是 V 代表你的 x,y

   V = atan2( y-y0 , x-x0 ) - ang0;