Orbital Mechanics, Vector Equation

hi everyone,

my goal is to make a simulation of the solarsystem in which i can fly around with a spaceship. the first problem is that the distances between object and the relation to their sizes are far too big for “usual” float calculations. to make the furthest body (lets say neptune) visible in my view frustum (zfar ~= 1000.0f), i’ll draw much smaller spheres much closer to me (the scene’s “camera”).

the problem is the mathematical equation to calculate a planets position for each time, like:
vec3 position = s(t);

a simlpe circle is the easiest thing, there wouldnt be changes in velocity over time, but an elliptic orbit is much more difficult. i’ve seached a few hours, could find anything useful (any formular like s(t), i found only r(phi) and v(phi)). can anyone give me a hint ?

the problem is that the velocity (and thus the angular velocity) isnt const, it varies depending on the distance to the center / sun. i’ve heared some time ago that just going with the force, calculating acceleration, and moving the object:

velocity += (Force / Mass) * frametime
position += velocity * frametime

… wont work because that makes the whole systen unstable, planets wont stay in an elliptic orbit.

any suggestions ?

There isn’t a closed-form expression for position as a function of time. Specifically, finding the eccentric anomaly E from the mean anomaly M requires solving Kepler’s equation, M=E-ε*sin(E), numerically, e.g. via Newton’s method.


double M = 2 * M_PI * (time % period) / period;
double E = M;
for (int i = 0; i < 5; i++)
    E -= (E - epsilon * sin(E) - M) / (1 - epsilon * cos(E))
double x = r_major * cos(E);
double y = r_minor * sin(E);

Having found E, the velocity can be found via the chain rule:

M=E-εsin(E)
=> dM/dt = dE/dt-ε
cos(E)dE/dt
= dE/dt * (1-ε
cos(E))
=> dE/dt = (dM/dt)/(1-ε*cos(E))
dM/dt = 2π/period
dx/dt = -sin(E) * r_major * (dE/dt)
dy/dt = cos(E) * r_minor * (dE/dt)

[QUOTE=GClements;1286462]There isn’t a closed-form expression for position as a function of time. Specifically, finding the eccentric anomaly E from the mean anomaly M requires solving Kepler’s equation, M=E-ε*sin(E), numerically, e.g. via Newton’s method.


double M = 2 * M_PI * (time % period) / period;
double E = M;
for (int i = 0; i < 5; i++)
    E -= (E - epsilon * sin(E) - M) / (1 - epsilon * cos(E))
double x = r_major * cos(E);
double y = r_minor * sin(E);

Having found E, the velocity can be found via the chain rule:

M=E-εsin(E)
=> dM/dt = dE/dt-ε
cos(E)dE/dt
= dE/dt * (1-ε
cos(E))
=> dE/dt = (dM/dt)/(1-ε*cos(E))
dM/dt = 2π/period
dx/dt = -sin(E) * r_major * (dE/dt)
dy/dt = cos(E) * r_minor * (dE/dt)[/QUOTE]

Wow! That is amazing! A lot of it goes over my head, especially the Calculus, but still it’s quite something even if I don’t have enough knowledge to fully appreciate it. I’m amazed that someone would know that. Where did you learn this?

Kepler’s equation can be found e.g. on Wikipedia. The rest is just basic calculus.

My suggestion is to get an orbit mechanics software library that you can link to.
They are available for free. I’m guessing they would have planetary propagators built in.
Google something like ‘Orbit Mechanics Software’, or ‘Planet Propagation Software’.
You want something simple that you can link in with your graphics program.
This way you can concentrate on the graphics without worrying about coding up propagation routines correctly.
Here’s an example -

     [https://sourceforge.net/projects/libnova/?source=directory](https://sourceforge.net/projects/libnova/?source=directory)

Good luck.

undestood! thank you all for your replies

@ Carmine: i’ve taken alook into that, interesting, but its faster to just make a function that does what i want instead figuring out how that lib works :smiley:

i didnt know that there isnt a solution to the equation, just to reiterate:
0. assumption is, that the positive X-axis points to the short end of the ellipse, where the velocity is the highest, right ? (to the “perihelion”, in the picture left, correct ? origin = sun position, correct ?)

  1. calculate M, which is just the mean “angular frequency”, using T which is given for each planet

  2. with M as the start value, solve the kepler equation with newton approximation:

M = E - e * sin(E)
<=>
0 = E - e * sin(E) - M

now, these are all constants right ? E, e and M dont depend on time …

then you’ve done 5 iteration to get an approximation of E, so first we have to build d(E - e * sin(E) - M) / dE to do that:

d(E - e * sin(E) - M) / dE = 1 - e * cos(E)

then newton’s iteration:
Enext = Eprevious - (Eprevious - e * sin(Eprevious) - M) / (1 - e * cos(Eprevious))

finally we can get the 2D position:
x = a * (cos(E) - e)
y = b * sin(E)

next assumption / question:
a = the longest “radius” in x-direction, and b = the “width” or radius in y-direction, is that correct ???
if so, then we have to use the r(phi) function first to get b, because for each planet we only have perihelion / aphelion …

https://en.wikipedia.org/wiki/Orbit#Newtonian_analysis_of_orbital_motion
r(phi) = a(1 - e²) / (1 - e * cos(phi))

r(phi = 90°) = b

… where a in this case is the middle of both, perihelion and aphelion, called “semi-major axis”

The point where E=M=0, x=r_major, y=0, corresponds to the periapsis (perihelion for orbit about the sun, perigee for orbit around the earth), i.e. the point of closest approach to the attractor, where speed is greatest.

M is the mean anomaly. For a circular orbit, this is just the angle between the T=0 reference direction and the moving body. In the general case, it’s just the orbital proportion (time divided by orbital period) converted to an angle. dM/dt is constant, i.e. M increases at a constant rate.

e (epsilon) is the eccentricity, which is constant. It’s equal to sqrt(1-b2/a2), where a is the semi-major axis and b is the semi-minor axis.

M and E are both angles which change with time. M changes at a constant rate. E is calculated from M, and is the angle in the parametric equation for an ellipse: x=acos(E), y=bsin(E).

So you calculate M as 2pi(t-t0)/T where t is the current time, t0 is the time at periapsis, T is the orbital period. Then calculate E from M by solving Kepler’s equation via Newton’s method (or bisection, or whatever).

Right; that gives you x relative to the focus (rather than the centre).

Yes. a is the semi-major axis, b is the semi-minor axis.

[QUOTE=john_connor;1286483]if so, then we have to use the r(phi) function first to get b, because for each planet we only have perihelion / aphelion …
[/QUOTE]
The attractor is at the focus of the ellipse. The distance of the focus from the centre is f=ae. The distance of the focus from the periapsis is a-f = a(1-e); the distance of the focus from the apoapsis is a+f=a*(1+e). The mean of those two is the semi-major axis, a. The ratio of their difference to their sum is the eccentricity, e. The semi-minor axis b can be found via b=a/sqrt(1-e2).

In short, a, b, e, f, a-f, a+f are all related; given any two, you can find the others. That much is just geometry.

An orbit requires two parameters to describe the ellipse, two to describe the orientation of the ellipse, and two more to define the period and phase. There are many different ways in which these (the orbital elements) can be specified.

OK, thx again, now i’m getting it …
it’s been a while since i’ve done math on an ellipse, and the physics part just confused me a lot :dejection: :smiley: (energy, force, mass, acceleration, velocity, blabla … the “key” here is the kepler equation M = E - esinE & basic ellipse geometry)

another question, same topic:

i can now calculate the 3D orbital position of a planet, but i’ve got some headache from trying to transform a “space” reference system into my “scene” system

imagine the folowing:

dvec3 PositionPlanetInSpace;
dvec3 PositionAstronaut;

struct {
vec3 Position;
quat Rotation;
float FieldOfView = radians(45.0f);
float Znear = 0.1f;
float Zfar = 1000.0f;
} Camera;

the center of the “space” coordinatesystem = dvec3(0, 0, 0). the “PositionPlanetInSpace” is calculated in that coordinatesystem, so the numbers can be very high (double). i want to look at that system from a certain (lets say static) position “PositionAstronaut”, that should represent the cameras position in “space”.

how would you solve the problem that NOTHING is visible ??
either the planet is too far away or clipped away because its not in the "Camera"s view frustum (Zfar = only 1000.0f)

should i just do all the shading with "double"s and increase the Zfar to … you now what …
i thought it would be easier to just multiply the dvec3 distance with a scalefactor, but … i dont sonehow get it right :frowning:

dvec3 asseenfromspace = planet.Position - PositionAstronaut;
double distance = length(asseenfromspace);
float scalingdown = Camera.Zfar * distance / MAX_DISTANCE;
/* MAX_DISTANCE = <double>farthest planet */
/* if only 1 planet there: distance / MAX_DISTANCE = 1 */
vec3 virtualposition = asseenfromspace * scalingdown;

that was my idea, however …

by the way, a sphere has radius = 1 in model space, and planets have “double Radius = 6371000.0”

OK, it seems that i forgot that the distances in the solar sytem is extremely huge, compared to the planet’s sizes:
for examle, if i’m sitting at the sun’s position and look to earth, the “view angle” is:

phi ~ tan(phi) = radius / distance = 6371km / 149598023km = 0.0000425rad = 0.002°

… which means that i have to scale them up, just to be able to see them (without a virtual telescope ^^)

Maybe this will also help you: http://www.stjarnhimlen.se/comp/ppcomp.html (sorry, the forum does not allow me to post links)