Нахождение истинной аномалии по векторам состояния - PullRequest
1 голос
/ 27 мая 2019

Я пытаюсь преобразовать векторы состояния (положение и скорость) в элементы Кеплера, однако я сталкиваюсь с проблемами, когда отрицательная скорость или положение приведут к неверным результатам при попытке вычислить истинную аномалию.

Вот несколько способов вычислить истинную аномалию:

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

Редактировать: еще один способ сделать это, на что указал Призрак:

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

Все позиции относятся к родительскому телу.

Все эти функции согласуются, если position.x, position.y и speed.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) };

Смотрите мой ответ для дополнительного кода о том, как я тестирую, и о том, как ям, используя для некоторых других переменных.

На этом я буду ходить кругами.это результат ошибки в моем существующем коде, которая обнаруживается при одних условиях, но не при других.Я предполагаю, что реальный вопрос сейчас - ПОЧЕМУ Получаю ли я другие результаты с позицией / скоростью выше, которые не соответствуют моим ожиданиям или друг другу?

Ответы [ 3 ]

1 голос
/ 27 мая 2019

Предполагая 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 Я вижу это так:

img

Поскольку вы получили свои координаты относительно родителя, вы можете предположить, что они уже находятся в центре фокальной точки, поэтому больше не нужно 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;
    
1 голос
/ 28 мая 2019

Хорошо, поэтому при дальнейшем тестировании кажется, что все мои исходные вызовы возвращают правильные значения, однако, когда я смотрел на выходные данные, я не принимал во внимание 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);
}

Спасибо Futurogogist и Spektre за помощь.

1 голос
/ 27 мая 2019

Я предполагаю, что вы работаете в двух измерениях?

Двумерные векторы положения 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);
}
...