Proyectar una elipse en un sistema de coordenadas en 3D

Buenas, a ver si me podéis echar un cable. Antes de nada os pongo en situación:

El año pasado en una asignatura de programación, en una de las prácticas hice un sistema de órbitas planetarias muy básico. Estaba programado en C++ usando solo canvas siendo por tanto un dibujo en dos dimensiones. Quedó bastante chulo pero era algo limitado, ya que por ejemplo, no podía tener en cuenta la inclinación orbital de los planetas.
Este año en una asignatura optativa llamada Informática Gráfica, seguimos programando en c++, solo que usamos APIS gráficas como OpenGL o DirectX para poder hacer cosillas como videojuegos o pequeñas simulaciones. Eso me da más libertad y puedo trabajar en un sistema de coordenadas de tres dimensiones mucho más fácilmente, por lo que me había propuesto reaprovechar la práctica del año pasado para practicar un poco, introduciendo dicho factor de inclinación orbital y ya veré si otros más mas adelante.

Hasta ahora he conseguido ya hacer un sistema de órbitas elípticas sobre un plano simple, por lo que decidí introducir el atributo de inclinación orbital. El principal problema es que aunque una órbita elíptica está sobre un plano, en el momento en el que inclinas ese plano, la elipse gana automáticamente una coordenada más que no sé como calcular, por lo que procedí a hacer un invento que funciona a medias:

Calculo primero la posición del planeta en el plano y a continuación inclino las coordenadas x grados para ponerlas en su posición (como si pones un punto con bolígrafo en la puerta, y a continuación la abres unos grados, el punto se desplaza con la puerta al nuevo plano que has abierto). El caso es que funciona solo a medias, cuando va por la mitad de la órbita (es decir, cuando el planeta llega a la posición donde el plano inclinado se corta con el plano no inclinado),  empiezan a suceder cosas raras.

¿Qué busco?
Podéis tratar de ayudarme de varias maneras:

1) Voy a pastear más abajo el código de la parte del programa que recalcula las posiciones del planeta, si tenéis alguna idea de algo que debería cambiar para que funcione y funciona eso, entonces premio para el que consiga encontrar el error que yo no veo.

2) Buscar una forma alternativa de calcular las coordenadas de un planeta en una órbita inclinada. Lo idóneo sería algo que trabaje con fórmulas matemáticas lo más parecidas posible a las ecuaciones paramétricas de la elipse pero que me permita trabajar con 3 coordenadas en vez de con 2, ya que actualmente esas fórmulas son las más cómodas para programar.

Código:
Código del programa

Algunas cosas a tener en cuenta:

- En mi sistema de coordenadas, la y es la vertical, la x la horizontal y la z la profundidad.
- El sol (que es el foco de la elipse) se encuentra situado en (x,y,z) = (0,0,0)
- Las fórmulas que uso son básicamente las paramétricas.
- Si tenéis alguna duda más del código, intentaré explicarlo más.
- Como me parece una cuestión bastante jodida, ofrezco una suculenta recompensa (500 gmfos) al que me dé una solución o descubra lo que falla en el código para que funcione. Si nadie da con una solución pero me dais ideas frescas para pensar en otras direcciones diferentes, repartiré gmfos entre las ideas que me sean más útiles para desarrollar un nuevo algoritmo.