I'm a software developer but my degree with in mathematics. I'm trying to develop a game/simulation of a fictitious randomly generated solar system, and trying to work out how to model and calculate the planetary orbits. I've been reading what resources I can find and understand, especially https://calgary.rasc.ca/orbits.htm and feel I am close to what I need, but stuck on a few points at the end.

Firstly I think I can happily ignore:

  • That the planet(s) orbit and the sun itself orbit around the solar system barycenter / that the gravity of the planet(s) pull on the sun. I can assume the sun is a fixed point (0,0,0) and the planets orbit around that.
  • That the gravity of the planet(s) influence each other.

Given that, what I plan to do for a given random planet so far is:

  1. Roll random value for "a", semi-major axis of elliptical orbit (in metres, for benefit of calculation below)
  2. Roll random value for "e", eccentricity, between 0 <= e <= 0.8 seems reasonable
  3. Calculate "b", semi-minor axis = a * √(1 - e2)
  4. Foci are at +/- ae so can now translate the ellipse ae to the right to place one of the foci at (0,0,0) where the sun is
  5. This is so far assuming the ellipse axes are lined up with the coordinate system axes so the ellipse is oriented horizontally, so now rotate the ellipse in the 2D plane a random number of degrees 0 <= r < 360 (I'm afraid I have not yet found what the correct terminology is for this)
  6. Roll random value for orbit inclination, between -20 <= i <= +20 seems reasonable, and rotate the ellipse in 3D space accordingly
  7. Calculate orbit duration T in seconds = 2π × √(a3 / GM) where M is the mass of the sun

Points on the ellipse can be calculated as x = a cos t, y = b sin t (then applying the translate/rotate operations in 4-6 above). So this is enough to draw the elliptical orbit. But clearly the points generated here don't take into account that the orbit is faster when the planet is closest to the sun vs its slowest when the planet is at its furthest point from the sun.

I can roll a random point on the ellipse as the starting point for the planet. But how can I calculate what point on the ellipse the planet will have each day? I see two possible approaches, I'm happy for either:

  • As a definitive calculation, so a formula into which I could plug in the parameters of the ellipse and calculate coordinates for any time "t". or
  • If I knew a starting point and velocity, I could simply advance the planet by that amount, and apply the effect of the sun's gravity on its velocity vector. (I wrote a simplistic demo that did just that many years ago). But this requires knowing the planet's position and direction of travel on day 1 (which I can work out as above) but also its speed to make it follow the prescribed orbit (which I do not know how to calculate, other than by a loop of guess, calculate, refine guess until I reach an acceptable level of accuracy).

Many thanks for any suggestions.

  • 1
    $\begingroup$ Kepler's second law says equal areas of the ellipse in equal amounts of time. Mean Anomaly and Eccentric Anomaly on the same page are what you're looking for as far as locations along an ellipse. You should do some more research on values for orbital parameters; for example, notably-eccentric Pluto's e is only 0.2488. $\endgroup$
    – Erin Anne
    Commented Apr 5, 2023 at 3:45

1 Answer 1


Wow, I think you are almost there!

Here are some goodies to get you up and running as a creator of planetary systems

  1. Roll random value for "e", eccentricity, between 0 <= e <= 0.8 seems reasonable

As @ErinAnne points out high eccentricities like 0.8 are pretty dramatic and should be used with a really really low probability, otherwise your solar system will be crazy-looking and very unphysical. On astronomical timescales high eccentricity objects either get ejected or (more frequently?) get "calmed down" due to those weak n-body interactions with other planets.

Roughly speaking, planetary systems can, under some conditions, "cool off" over time.

You haven't mentioned masses (because you'll use only Keplerian orbits and won't do n-body) but in the thousand(s?) of known exo-planetary systems catalogued there are usually one or more often a few Jupiter-sized objects. These gravitational "moderators" (or perhaps "bullies") are constantly telling folks to "Settle down people! Find your place, take a seat" and often help planets tend towards some level of circularity.

So rather than a flat random distribution of eccentricity between 0 and "really high", I recommend either a normal (Gaussian) distribution centered at 0 or say 0.1, or a pseudo-thermal (exponentially decreasing) distribution of eccentricities.

If you'd like some suggestions how to use random number generator libraries to get truncated normal (Gaussian) distributions or thermal distributions (or others) please let me know and I'll add some here.

Or you can also ask that as a new question. You might consider Scientific Computing SE as long as perhaps you call it a "simulation" and not a "game". ;-)

I can roll a random point on the ellipse as the starting point for the planet. But how can I calculate what point on the ellipse the planet will have each day?

You could do numerical integration but I guess you would want to avoid that. If you don't want to avoid that then let me know and I'll add links to several helpful posts. For completeness I'll add a short script at the end. It uses a standard library integrator but you can "roll your own". There are some flavors of Euler that are fast and accurate, or a variable step size RK45.

But the biggest problem with numerical integration is that it always gets worse over time!** So the series expansion or inverse solution + interpolation (see below) are preferable because they are strictly periodic - you always subtract an integer number of periods $T$ before calculation the position within one orbit!

Instead, what's really cool to find out is that for a Keplerian orbit (elliptical or otherwise) while position as a function of time can't be calculated analytically (you need a series approximation or an integration) the inverse problem: time as a function of position - has an analytical solution! See for example:

So you can either use a series approximation for the forward problem (position as a function of time) which will converge quickly for modest eccentricities (another reason to avoid those "crazy point-eights") and need more terms for the high eccentricities, or choose and equal or unequally spaced set of angles (called eccentric anomalies in Kepler's space words) use it to get $\mathbf{r}$ and time, then interpolate that with your own daily timepoints.

(Feel free to ask a new question about the series approximations, but I think it's discussed in one Wikipedia article or another. I just can't find it right now.)

For orbital inclination and longitude of ascending node

These two Keplerian orbital parameters define the plane of the orbit. The best way to generate pairs randomly in my opinion (or at least the coolest) is to appreciate that each orbital plane has a single normal vector, and they all point from your orbital center (your Sun in your scheme) to some place on a unit sphere.

Retrograde orbits (circling the "wrong way" are defined by points in the lower hemisphere of the same unit sphere.

Get out your favorite "random points on a sphere" algorithm to choose these points. Their $\theta$ (similar to latitude but measured from the top of the sphere in spherical coordinates) will be the orbit's inclination, and their $\phi$ (similar to longitude) will be the orbit's longitude of ascending node.

If you don't already have a favorite, see the algorithms, links and code in all of the answers to

Mathematically related and potentially interesting to Math majors dabbling in orbital mechanics for now, hopefully you'll dive in deeper!:

Just in case you want to get the positions by numerical integration, here's a short, simple 2D Python script to get you started. It's not meant as a solution, it just gives a "taste" of how a numerical solution would work.

numerical integration of a 2D orbit with reduced and physical units (click for full size)

Python script:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import time

def deriv(t, y, GM): # time derivative of state vector
    x, v = y.reshape(2, -1)
    acc = -GM * x * ((x**2).sum())**-1.5  # - GM * x_hat / abs(x)^2
    return np.hstack((v, acc))

two_pi = 2 * np.pi

e = 0.8   # eccentricity

# reduced units
GM_red = 1.   # standard gravitational parameter 
a_red = 1.   # semimajor axis

# physical units
GM_phys = 1.3271244001E+20 # m^3/s^2 (Sun) 
a_phys = 149598023000. # meters (for Earth's orbit)

# okay go for it!
names = 'reduced units', 'physical units'
semis = a_red, a_phys
GMs = GM_red, GM_phys

method = 'DOP853'
dense_output = False
n_eval = 201
tol = 1E-08

answers = []
for (a, GM, name) in zip(semis, GMs, names):
    # initialize state vector y0
    r_peri = (1 - e) * a
    v_peri = (GM * (2/r_peri - 1/a))**0.5  # https://en.wikipedia.org/wiki/Vis-viva_equation
    y0 = np.array([r_peri, 0, 0, v_peri], dtype=float)

    # initialize evaluation grid & spans (interpolator after integration finishes) 
    T_period = two_pi * (a**3/GM)**0.5 # https://en.wikipedia.org/wiki/Elliptic_orbit
    t_eval = np.linspace(0, T_period, n_eval)
    t_span = t_eval.min(), t_eval.max()

    args = (GM, )
    print('starting: ', name)
    t_start = time.process_time()
    print('t_span: ', t_span)
    print('y0: ', y0)
    print('args: ', args)

    answer = solve_ivp(deriv, t_span, y0, method=method, t_eval=t_eval,
                       dense_output=dense_output, events=None, args=args, 
    print('finished integrating, process time (sec): ',
          time.process_time() - t_start)

# okay plot the orbits
if True:
    fig, axes = plt.subplots(2, 1, figsize=[10, 6])
    ax1, ax2 = axes
    titles = 'reduced units', 'physical units'
    for answer, ax, name, title in zip(answers, axes, names, titles):
        x, y, vx, vy = answer.y
        ax.plot(x, y, '-k')
        ax.scatter(x[::5], y[::5], marker='o', c=t_eval[::5], cmap='jet')
        ax.plot([0], [0], 'oy', ms=20)
        difference = answer.y[:, -1] - answer.y[:, 0]
    plt.suptitle('eccentricity: ' + str(e))
  • 2
    $\begingroup$ inclination to the ecliptic is also pretty tightly-bounded for large bodies (i.e. planets) in our solar system. A planetary system with orbital planes chosen by the random-vector-on-a-sphere is going to have some flavor to it. I don't know if it's contrary to physics though (I thought accretion discs tended to cancel out out-of-ecliptic velocity but that's definitely not my wheelhouse) $\endgroup$
    – Erin Anne
    Commented Apr 5, 2023 at 6:07
  • $\begingroup$ @ErinAnne That's a good point! Inclinations should also probably be generated from some normal distribution as well. While I'm probably all "exchanged out" for the day I will try to add that in the next edit. You're more than welcome to beat me to it and make an edit yourself. $\endgroup$
    – uhoh
    Commented Apr 5, 2023 at 7:21
  • 1
    $\begingroup$ Many thanks for the detailed answers. I should have plenty to go from there. My range of values for eccentricity was from looking at diagrams like on the page Erin linked - if you look at the ellipse with e=0.2 it does not obviously strike you as an ellipse, all you can see is that the sun is off-centre. If they are so close to all just being circles then I may as well not bother and just make all orbits circular. But agree 0.8 is probably way too high as it places the planet far too close to the sun at its closest point. $\endgroup$
    – nigelg
    Commented Apr 5, 2023 at 10:16
  • $\begingroup$ @nigelg aesthetics is important for your application, the reason Mars' small 0.09 eccentricity is as obvious as it is, is that we compare it to nearby Earth's nearly circular orbit. So your generating algorithm might use some simple rules to throw out "boring" or "crazy-looking" systems. I don't think you need rigorous randomness as much as you need your systems to have a bit of "oh that looks cool" personality. Go make a thousand and then start deciding which cases you'd like to avoid, and have fun either "playing god" or Worldbuilding :-) $\endgroup$
    – uhoh
    Commented Apr 5, 2023 at 10:28
  • 1
    $\begingroup$ @nigelg Sure, at e=0.2, it's hard to see that the orbit isn't a circle, apart from the displacement of the Sun. But there is a noticeable affect on the orbit speed. In this plot, the 36 dots are placed at equal time steps. i.sstatic.net/g2saC.png $\endgroup$
    – PM 2Ring
    Commented Apr 5, 2023 at 23:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.