1

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

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

/// <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) };

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

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

4

3 回答 3

2

假设 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
    

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

请参阅[Edit3] 进一步提高 Newton D'ALembert 积分精度

有关更多信息和方程式,请参见:

[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. 初始ang0角度a

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

    ang0 = atan2( ya-cy , xa-cx );
    
  4. 最后V是你们中的任何一个x,y

    V = atan2( y-y0 , x-x0 ) - ang0;
    
于 2019-05-27T07:04:48.453 回答
2

好吧,进一步测试看来,我的原始计算确实都返回了正确的值,但是当我查看输出时,我没有考虑 LoP 并且基本上没有认识到 180 与 -180 的角度基本相同。(我也在查看弧度的输出,只是没有看到应该很明显的内容) 长话短说,我有一个错误,我认为是在代码的这个区域并且迷失在杂草中。
看来我上面错了。有关边缘情况,请参见 OP。
这是我用来测试这些的一些代码,
我使用了以下输入的变体:

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

并用以下方法测试它们

double parentMass = 1.989e30;
double objMass = 2.2e+15;
double sgp = GameConstants.Science.GravitationalConstant * (parentMass + objMass) / 3.347928976e33;

Vector4 ev = OrbitMath.EccentricityVector(sgp, pos, vel);
double e = ev.Length();
double specificOrbitalEnergy = Math.Pow(vel.Length(), 2) * 0.5 - sgp / pos.Length();
double a = -sgp / (2 * specificOrbitalEnergy);
double ae = e * a;
double aop = Math.Atan2(ev.Y, ev.X);
double eccentricAnomaly = OrbitMath.GetEccentricAnomalyFromStateVectors(pos, a, ae, aop);
double aopD = Angle.ToDegrees(aop);
double directAngle = Math.Atan2(pos.Y, pos.X);

var θ1 = OrbitMath.TrueAnomaly(sgp, pos, vel);
var θ2 = OrbitMath.TrueAnomaly(ev, pos, vel);
var θ3 = OrbitMath.TrueAnomalyFromEccentricAnomaly(e, eccentricAnomaly);
var θ4 = OrbitMath.TrueAnomalyFromEccentricAnomaly2(e, eccentricAnomaly);
var θ5 = OrbitMath.TrueAnomaly(pos, aop);
double angleΔ = 0.0000001; //this is the "acceptable" amount of error, really only the TrueAnomalyFromEccentricAnomaly() calcs needed this. 
Assert.AreEqual(0, Angle.DifferenceBetweenRadians(directAngle, aop - θ1), angleΔ);
Assert.AreEqual(0, Angle.DifferenceBetweenRadians(directAngle, aop - θ2), angleΔ);
Assert.AreEqual(0, Angle.DifferenceBetweenRadians(directAngle, aop - θ3), angleΔ);
Assert.AreEqual(0, Angle.DifferenceBetweenRadians(directAngle, aop - θ4), angleΔ);
Assert.AreEqual(0, Angle.DifferenceBetweenRadians(directAngle, aop - θ5), angleΔ);

和以下比较角度:

public static double DifferenceBetweenRadians(double a1, double a2)
{
    return Math.PI - Math.Abs(Math.Abs(a1 - a2) - Math.PI);
}

并因此找到偏心率向量:

public static Vector4 EccentricityVector(double sgp, Vector4 position, Vector4 velocity)
{
    Vector4 angularMomentum = Vector4.Cross(position, velocity);
    Vector4 foo1 = Vector4.Cross(velocity, angularMomentum) / sgp;
    var foo2 = position / position.Length();
    return foo1 - foo2;
}

和偏心异常:

public static double GetEccentricAnomalyFromStateVectors(Vector4 position, double a, double linierEccentricity, double aop)
{
    var x = (position.X * Math.Cos(-aop)) - (position.Y * Math.Sin(-aop));
    x = linierEccentricity + x;
    double foo = GMath.Clamp(x / a, -1, 1); //because sometimes we were getting a floating point error that resulted in numbers infinatly smaller than -1
    return Math.Acos(foo);
}

感谢 Futurogist 和 Spektre 的帮助。

于 2019-05-28T09:19:58.213 回答
1

我假设你在二维工作?

p位置和速度的二维向量v。该常数K是重力常数与重力产生体质量的乘积。计算偏心率向量

eccVector = (dot(v, v)*p - dot(v, p)*v) / K - p / sqrt(dot(p, p));

eccentricity = sqrt(dot(eccVector, eccVector));

eccVector = eccVector / eccentricity;

b = { - eccVector.y, eccVector.x};  //unit vector perpendicular to eccVector  

r = sqrt(dot(p, p));

cos_TA = dot(p, eccVector) / r;   \\ cosine of true anomaly
sin_TA = dot(p, b) / r;           \\ sine of true anomaly

if (sin_TA >= 0) { 
   trueAnomaly = arccos(cos_TA);
}
else if (sin_TA < 0){
   trueAnomaly = 2*pi - arccos(cos_TA);
}
于 2019-05-27T02:07:30.457 回答