Actualizado 2016-12-12

Estoy escribiendo un programa de mecánica orbital en Unity . Convierto la posición y la velocidad de un objeto en órbita en elementos orbitales, luego convierto los elementos orbitales nuevamente en vectores de posición cartesianos para poder trazar la órbita completa. Seguí las ecuaciones en estos dos enlaces:

Vectores de estado cartesiano a elementos keplerianos

Elementos keplerianos a vectores de estado cartesianos

Las ecuaciones que estoy siguiendo para calcular los vectores de estado suponen que el plano de referencia es XY. En Unity, el plano de referencia es XZ, siendo Y la "dirección hacia arriba". Entonces, lo que hice fue restar 90 grados de la inclinación, de modo que una órbita en el plano XZ tendría una inclinación cero, lo que terminó siendo la causa de mis problemas. Así que eliminé esa modificación y ahora la mayoría de mis órbitas se muestran correctamente, para las órbitas que comienzan en el periapsis. Pero las órbitas que comienzan en la apoapsis están invertidas. Vea las imágenes a continuación:

Aquí hay una órbita correctamente dibujada. Tengo una pequeña nave espacial en órbita y comienza en el eje -X. Dibujé una órbita incompleta para poder ver el punto de inicio de la órbita, que puede ver que comienza en el eje -X.

Aquí cambié el estado inicial de la nave para que comience en apoapsis. Como puede ver, el punto de partida de la órbita está a 180 grados, y se encuentra en el eje +X en lugar del eje -X. De hecho, la órbita siempre se dibuja a partir del periápside. Lo que estoy tratando de hacer es que el primer punto de la órbita sea la posición de la nave en el momento en que llamo a la función ConvertToCartesian

Aquí está mi código:

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

Para dibujar la línea de la órbita, simplemente llamo a la función ConvertToCartesian en un ciclo con valores crecientes de t para obtener una matriz de vectores de posición.

Valores: M = 3.3e16, GM ~= 2.2e6, r = (-50000, 0,0) Para una órbita circular, v = (0,0, sqrt(GM/50000))

Pregunta de seguimiento: dado que las ecuaciones asumen el plano de referencia XY, ¿cómo puedo cambiarlo para que se ajuste al plano de referencia XZ que estoy usando en mi programa?

Actualizar

Actualicé el método ConvertToCartesian, de modo que toma la anomalía verdadera (nu) calculada a partir de la función ConvertToKeplerElements como entrada. Luego calcula la anomalía excéntrica (E), luego calcula la anomalía media (Mt) a partir de eso. Cuando se llama a ConvertToCartesian en bucles posteriores, Mt se actualiza con el paso de tiempo relevante. Luego calculo E nuevamente a partir del Mt actualizado, luego calculo nu a partir de la E actualizada. Esto resuelve la mayoría de los casos. El problema es que, a veces, nu se calcula a 180 grados de donde debería estar, lo que hace que el punto de inicio de la órbita mostrada esté a 180 grados de distancia.

En mi función ConvertToKeplerElements, de acuerdo con las matemáticas, nu = 2pi - nu si rv < 0. Esto sucede, por ejemplo, cuando configuro los vectores r, v iniciales de modo que la nave no comience ni en el periapsis ni en el apoapsis, sino en algún punto intermedio. Sé que esto es correcto, pero cuando calculo nu->E->nu, ese giro de 180 grados se pierde en el cálculo.

Aquí están las partes relevantes del método actualizado 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)); 

Aquí está el mal resultado que estoy obteniendo:

La verdadera anomalía al comienzo de la órbita debería ser 1,5pi (~4,71), pero al calcularla desde E es 0,5pi (~1,57)