Обновлено 2016-12-12

Я пишу программу орбитальной механики в Unity . Я преобразовываю положение и скорость объекта на орбите в элементы орбиты, а затем преобразую элементы орбиты обратно в декартовы векторы положения, чтобы я мог построить всю орбиту. Я следил за уравнениями по этим двум ссылкам:

Декартовы векторы состояния в кеплеровы элементы

Кеплеровы элементы в декартовых векторах состояния

Уравнения, которым я следую для расчета векторов состояния, предполагают, что эталонная плоскость — это XY. В Unity эталонной плоскостью является XZ, а Y — «направление вверх». Итак, что я сделал, так это вычел 90 градусов из наклонения, чтобы орбита в плоскости XZ была с нулевым наклонением, что в итоге и стало причиной моих проблем. Поэтому я удалил эту модификацию, и теперь большинство моих орбит отображаются правильно для орбит, начинающихся в перицентре. Но орбиты, начинающиеся в апоапсисе, обратные. Смотрите изображения ниже:

Вот правильно нарисованная орбита. У меня есть небольшой космический корабль на орбите, и он начинается с оси -X. Я нарисовал неполную орбиту, чтобы увидеть начальную точку орбиты, которая, как вы видите, начинается на оси -X.

Здесь я изменил начальное состояние корабля, чтобы он стартовал в апоапсисе. Как видите, начальная точка орбиты смещена на 180 градусов и лежит на оси +X, а не на оси -X. На самом деле орбита всегда рисуется, начиная с перицентра. Я пытаюсь сделать так, чтобы первая точка орбиты была позицией корабля в момент времени, когда я вызываю функцию ConvertToCartesian.

Вот мой код:

Vector3 ConvertToCartesian(float e, float a, float i, float O, float w, float t, float t0){
    float mu = G * M;
    float Mt; //Mean anomaly at time t
    float T = 2 * Mathf.PI * Mathf.Sqrt(a * a * a / mu); //Orbital period
    if (t == t0) {
        t = t0;
        Mt = 0; 
    } else {
        float deltaT = (T/60) * (t - t0); //divide time increments into 1/60th of the orbital period
        Mt = deltaT * Mathf.Sqrt (mu / Mathf.Pow (a, 3));
    }
    //Calculate eccentric anomaly using Newton's method
    float E = Mt;
    float F = E - e * Mathf.Sin (E) - Mt;
    int j = 0, maxIter = 30;
    float delta = 0.000001f;
    while (Mathf.Abs(F) > delta && j < maxIter) {
        E = E - F / (1 - e * Mathf.Cos (E));
        F = E - e * Mathf.Sin (E) - Mt;
        j++;
    }
    float nu = 2 * Mathf.Atan2 (Mathf.Sqrt (1 + e) * Mathf.Sin (E / 2), Mathf.Sqrt (1 - e) * Mathf.Cos (E / 2)); //True anomaly

    float rc = a * (1 - e * Mathf.Cos(E)); //Distance to central body

    Vector3 o = new Vector3(rc * Mathf.Cos(nu), rc * Mathf.Sin(nu), 0); 
                
    Vector3 odot = new Vector3 (Mathf.Sin (E), Mathf.Sqrt (1 - e * e) * Mathf.Cos (E), 0); //Velocity vector in the orbital frame
    odot = (Mathf.Sqrt (mu * a) / rc) * odot;

    float rx, ry, rz;
    rx = o.x; ry = o.y; rz = o.z;

    rx = ( o.x * (Mathf.Cos (w) * Mathf.Cos (O) - Mathf.Sin (w) * Mathf.Cos (i) * Mathf.Sin (O)) -
            o.y * (Mathf.Sin (w) * Mathf.Cos (O) + Mathf.Cos (w) * Mathf.Cos (i) * Mathf.Sin (O)));
    ry = (o.x * (Mathf.Cos (w) * Mathf.Sin (O) + Mathf.Sin (w) * Mathf.Cos (i) * Mathf.Cos (O)) +
        o.y * (Mathf.Cos (w) * Mathf.Cos (i) * Mathf.Cos (O) - Mathf.Sin (w) * Mathf.Sin (O)));
    rz = (o.x * (Mathf.Sin (w) * Mathf.Sin (i)) + o.y * (Mathf.Cos (w) * Mathf.Sin (i)));

    Vector3 r = new Vector3(rx, ry, rz); //Position vector

    return r / objectScale;

}


/* Convert a body's cartesian state vectors (position, r, and velocity, v) into Kepler elements
 */
void ConvertToKeplerElements(Vector3 r, Vector3 v){
    float mu = G * M;
    Vector3 h = Vector3.Cross (r, v); //Orbital momentum
    Vector3 n = Vector3.Cross (Vector3.forward, h); 
    eVector = ((v.magnitude * v.magnitude - mu / r.magnitude) * r - Vector3.Dot(r, v) * v) / mu; //Eccentricity
    float emag = eVector.magnitude;
    float E = v.magnitude * v.magnitude / 2 - mu / r.magnitude; //Specific mechanical energy
    if (eVector.magnitude != 1) {
        a = -mu / (2 * E); //Semi major axis
        float p = a * (1 - emag * emag);
    } else {
        //a = infinity
        a = 0;
        float p = h.magnitude * h.magnitude / mu;
    }

    inc = Mathf.Acos (h.z / h.magnitude); //Inclination

    //Longitude of ascending node (LAN), angle of body to node vector
    if (inc == 0 || inc == Mathf.PI) { 
        O = 0; //For an equatorial orbit (i = 0), LAN is undefined. By convention, set to 0. 
    } else {
        O = Mathf.Acos (n.x / n.magnitude); 
    }
    if (n.y < 0) {
        O = 2 * Mathf.PI - O;
    }

    //Argument of periapsis
    if (eVector.magnitude == 0) {
        w = 0; //For a circular orbit, by convention place w at ascending node (w = 0)
    } else {
        if (inc == 0 || inc == Mathf.PI) {
            w = Mathf.Atan2 (eVector.y, eVector.x);
        } else {
            w = Mathf.Acos (Vector3.Dot (n, eVector) / (n.magnitude * eVector.magnitude)); 
        }
    }
    if (eVector.z < 0) {
        w = 2 * Mathf.PI - w;
    }

    //True anomaly
    if (eVector.magnitude == 0) {
        if (inc == 0 || inc == Mathf.PI) {
            nu = Mathf.Acos (r.z/ r.magnitude); //For a circular non-inclined orbit, nu is angle of body from the x axis
            //(True Longitude)
        } else {
            //For a circular inclined orbit, nu is the angle between ascending node position vectors, (Argument of Latitude)
            nu = Mathf.Acos (Vector3.Dot (n, r) / (n.magnitude * r.magnitude));
        }
    } else {
        nu = Mathf.Acos (Vector3.Dot (eVector, r) / (eVector.magnitude * r.magnitude)); //True anomaly is angle between eccentricity and position vectors
        if (Vector3.Dot (r, v) < 0) {
            nu = 2 * Mathf.PI - nu;
        }
    }
}

Чтобы нарисовать линию орбиты, я просто вызываю функцию ConvertToCartesian в цикле с увеличением значения t, чтобы получить массив векторов положения.

Значения: M = 3,3e16, GM ~= 2,2e6, r = (-50000, 0,0) Для круговой орбиты v = (0,0, sqrt(GM/50000))

Дополнительный вопрос: поскольку уравнения предполагают базовую плоскость XY, как я могу изменить ее, чтобы она соответствовала базовой плоскости XZ, которую я использую в своей программе?

Обновлять

Я обновил метод ConvertToCartesian, чтобы он принимал в качестве входных данных истинную аномалию (nu), вычисленную функцией ConvertToKeplerElements. Затем он вычисляет эксцентрическую аномалию (E), а затем вычисляет среднюю аномалию (Mt). Когда ConvertToCartesian вызывается в последующих циклах, Mt обновляется с соответствующим временным шагом. Затем я снова вычисляю E из обновленного Mt, затем вычисляю nu из обновленного E. Это решает большинство случаев. Проблема в том, что иногда nu вычисляется на 180 градусов от того места, где оно должно быть, что делает начальную точку отображаемой орбиты отклоненной на 180 градусов.

В моей функции ConvertToKeplerElements, согласно математике, nu = 2pi - nu, если rv < 0. Это происходит, например, когда я задаю начальные векторы r, v так, что корабль стартует не в перицентре и не в апоапсисе, а где-то между ними. Я знаю, что это правильно, но когда я вычисляю nu->E->nu, этот переворот на 180 градусов теряется при вычислении.

Вот соответствующие части обновленного метода ConvertToCartesian.

        float mu = G * M;
    float E = Mathf.Acos((e + Mathf.Cos(nu)) / (1 + e * Mathf.Cos(nu)));
    float Mt = E - e * Mathf.Sin (E);

    //Mean anomaly at time t
    float T = 2 * Mathf.PI * Mathf.Sqrt(a * a * a / mu); //Orbital period
    float deltaT = (T/60) * (t - t0); //divide time increments into 1/60th of the orbital period
    Mt = Mt + deltaT * Mathf.Sqrt (mu / Mathf.Pow (a, 3));

    //Calculate eccentric anomaly using Newton's method
    float F = E - e * Mathf.Sin (E) - Mt;
    int j = 0, maxIter = 30;
    float delta = 0.000001f;
    while (Mathf.Abs(F) > delta && j < maxIter) {
        E = E - F / (1 - e * Mathf.Cos (E));
        F = E - e * Mathf.Sin (E) - Mt;
        j++;
    }
        
    //True anomaly
    nu = 2 * Mathf.Atan2 (Mathf.Sqrt (1 + e) * Mathf.Sin (E / 2), Mathf.Sqrt (1 - e) * Mathf.Cos (E / 2)); 

Вот плохой результат, который я получаю:

Истинная аномалия в начале орбиты должна быть 1,5pi (~4,71), но при расчете по E она составляет 0,5pi (~1,57)