Estoy modelando órbitas planetarias usando ecuaciones paramétricas, con el Sol en (0,0).
En 2D, es bastante sencillo ($a$ es la longitud del semieje mayor y $e$ es la excentricidad)
$$ x\left(t\right)=a\ cos\left(t\right)+ae,\\ y\left(t\right)=a\sqrt{1-e^2}\ sin\left(t\right)\\$$
Sin embargo, estoy teniendo dificultades para extender el modelo a 3D. Si utilizo coordenadas esféricas polares, $\phi$ es el ángulo entre el eje Z y el vector $\vec{OP}$, donde P es la posición del planeta.
Dado que la inclinación orbital se mide desde el eje x positivo hasta $\vec{OP}$, entonces $i=\frac{\pi}{2}-\phi$. De este diagrama, $z=\frac{r}{tan(\phi)}$.
Si utilizo la ecuación polar $r=\frac{a\left(1-e^2\right)}{1+ecos\ \theta}$, entonces $$z\left(t\right)=\frac{a\left(1-e^2\right)}{\left(1+ecost\right)tan\left(\frac{\pi}{2}-i\right)}=\frac{a\left(1-e^2\right)tan\left(i\right)}{\left(1+ecost\right)}$$
Sin embargo, al graficar la órbita de Venus y la Tierra en el mismo eje, obtengo este resultado en MATLAB. No creo que haya errores en mi código, por lo que debe ser un error en mis cálculos matemáticos.
tl;dr ¿Existe una mejor manera de modelar elipses con un ángulo específico de inclinación?
function [] = orbitEquations_3d(planetName, eccentricity, a, inclination)
% inicializar (1-e^2) y sqrt(1-e^2)
oneMinusE_sq = (1-eccentricity^2);
sqrtOneMinusE_sq = sqrt(oneMinusE_sq);
% convertir el ángulo de inclinación de grados a radianes
inclination = inclination/180 * pi;
% nombre del planeta
disp(planetName);
% mostrar ecuaciones
disp("paramétrico");
fprintf("x(t) = %.6f*cos(t) + %.6f\n", a, eccentricity*a);
fprintf("y(t) = %.6f*sin(t)\n", a * sqrtOneMinusE_sq);
fprintf("z(t) = %.6f*tan(%.6f)/(1+%.6fcos(t)) \n", a * oneMinusE_sq, inclination, eccentricity);
fprintf("\n\n");
end
