I was trying to fake elliptical orbits a few weeks ago. It turned out to be trivial to implement, but it took me forever to sort through lots of incomplete explanations and wrap my head around the math, so I thought I'd write it up...
We're going to specify the ellipse by giving its semi-major axis (half of the long axis) as a vector (
a = (ax, ay)) and the eccentricity e (0 is a circle, 1 is completely flattened, Mercury's orbit is the most eccentric in our solar system at e=0.2056, don't give e>0.9 to this algorithm). I give
a as a vector so that ellipses can be oriented in any direction.
We also need the orbital period
T (how long it takes for the planet to make one revolution around the star). Pick this to be whatever looks good to you. If you'd like the relationship between different orbits to be somewhat realistic, calculate
k=T*sqrt(star_mass/length(a)^3) for your sample orbit, and then use that
k in
T=k*sqrt(length(a)^3/star_mass) to calculate all your orbital periods.
We're going to start by pretending that the planet is moving around a circle at a constant speed. We'll call the angle to this pretend planet the
mean anomaly and calculate it as
M=2*PI*time/T (in radians, that is).
If you take the actual position of the planet and draw a line straight up to where it hits the
major circle, the angle from the center of the ellipse to
that point is the
eccentric anomaly E, as shown in this image from Wikipedia:
The mean anomaly and the eccentric anomaly are related by the equation
M = E - e*sinE. We have
M, so we want to solve this for
E. There's no closed-form solution, but it's simple to solve by Newton's method:
E = M // initial guess
do {
step = (E - e*sin(E)) / (1 - e*cos(E))
E += step
} while(abs(step) > small_number)
Now that we have
E, we can find the planet's position:
// star is at focus, we want ellipse center
c = star - e*a
// positions on "unit ellipse"?
x = cos(E) - e
y = sin(E) * sqrt(1 - e*e)
// scale and orient by our 'a' vector
a_perp = (-ay, ax) // or (ay, -ax) depending on which way your y-axis goes
planet = center + x*a + y*a_perp
For 3D I imagine you would just set
a_perp to be the cross product of
a with the normal vector for the plane of the ecliptic, or something like that.