Generating Orbital Elements from Abstract Solar Hierarchies

This post builds on the concepts contained in “Generating Random Simplex Solar Hierarchies” and the concepts contained in “The Keplerian Orbital Elements and Orbital Perturbations.”


In our quest to generate and visualize N-ary Simplex Solar systems, we have managed to generate an abstract simplex solar hierarchy, as shown in this previous post.

This abstract hierarchy contains no geometric data, only the binary tree of orbital nodes.

We wish to translate this abstract tree into a series of legitimate Orbital Elements, whose ellipses we want to render to screen.

Translating an Abstract Solar Hierarchy into Keplerian Elements

Recall from our previous post, “Generating Random Simplex Solar Hierarchies,” that our solar system hierarchies are comprised of a series of Orbital Node objects. Each Orbital Node has a reference to a parent node and two child nodes. This is reminiscent of the Quad Tree data structure we created for our Icospherical World Model, except that in this case, our data structure is a Binary Tree.

In all solar hierarchies, there will be an L0 (level 0) Orbital Node; the “root node.”

If this node has no children, then this is a unary solar system with one sun, and we can bypass this entire process and render that star system (which we will discuss in a later post).

However, if this node has children, then this node is a barycenter, and its two children might be either star systems, or barycenters themselves.

Determining L1 Separation

Starting at L1, the two children of the root node will need to have a distance which they are separated. Astronomically, the distance between two stars of a binary star system varies wildly – from less than one AU to more than a lightyear.

It turns out that this “separation radius” is normally distributed in log a, centered around a typical separation of tens of AU. If this means nothing to you, most binary star systems have a separation radius of 10 – 60 AU.

Typical Separation of Binary Star Systems; color corresponds to detection method

In the future, perhaps in PISES we can determine L1 separation in a more thoughtful, meaningful, astronomical way – but for now, just to get moving, let’s make L1 separation related to:

  1. How deep the solar hierarchy is (deeper systems will have a greater L1 separation)
  2. Some degree of randomness

In our SolarSystem constructor, we can determine the L1 separation like so:

// Generate Orbital Hierarchy (see previous post)
generateOrbitalHierarchy(); // this populates "allNodes"

// Obtain the Root Node
rootNode = allNodes.FirstOrDefault(n => n.depth == 0);

// Obtain the maximum depth of the hierarchy.
int maximumDepth = 0;
if (allNodes.Count > 1) 
    maximumDepth = allNodes.Aggregate((n1, n2) => n1.depth > n2.depth ? n1: n2).depth;

// Determine the L1 separation in AU.
L1Separation = maximumDepth * Numericals.randomGaussian(20, 6);

This is our first encounter with the Numericals class. In this class, I have a series of static mathematical helper methods. “RandomGaussian(mean, deviation)” does exactly what it sounds like it does: it takes in the arguments mean and deviation and returns a random value along the normal distribution centered on mean, with a standard deviation of deviation.

Determining LN Separation

Getting back to our code, this means that for binary star systems, the L1 separation will be (1 * randomGaussian(20, 6). For trinary systems, which must have a depth of 2 (with the root node at depth 0), our L1 separation will be (2 * randomGaussian(20, 6). The deeper the hierarchy, the greater the L1 separation.

The reason we are making deeper hierarchies physically larger is that the separation of L2, L3, and LN pairs is, astronomically, usually around [separation( L(N-1) ) / 2]. So, in the interest of not having our LN systems become squashed at their deeper levels, we make the L1 separation appropriately larger.

Now, with an L1 separation in hand, we can generate the orbits of the L1 sibling pair by calling the recursive method, “GenerateSiblingKeplerianOrbits().”

// Generate Keplerian Orbits for each node in the orbital hierarchy 
private void generateSiblingKeplerianOrbits(OrbitalNode node, int maximumDepth, float maximumSeparation)
    // If this was called on the root node, check for children. 
    if (node.depth == 0)
        // If there are children, recurse.
        if (node.child1 != null)
        {generateSiblingKeplerianOrbits(node.child1, maximumDepth, maximumSeparation);}

        // Otherwise return.

    // Otherwise, this is a node down in the hierarchy, and it must have a sibling. 
        // Obtain Sibling 
        OrbitalNode sibling = node.sibling;

        // Determine Separation Distance.
        float separation;

        // If this is an L1 pair, the separation distance is equal to the maximum separation.
        if (node.depth == 1)
        {separation = maximumSeparation;}
        // Otherwise, it's equal to half of the parent pair's separation.
        {separation = node.parent.barycenterSeparation / 2;}

    // ...

Accounting for Mass; Determining the Semimajor Axes

Now that we have obtained the separation distance for the sibling pair, we need to obtain the separation for each individual child of the sibling pair. This is because the children will be of differing masses; therefore, they will have different distances to their shared center of mass (their barycenter). These distances will eventually become the semimajor axes of the respective systems.

We can calculate the pair’s center of mass (in a one-dimensional scale, with ‘node’ at 0 and ‘sibling’ at maximumSeparation) and then determine the offset of both the node and its sibling, by continuing our “GenerateSiblingKeplerianOrbits()” method:

    // Determine the Offset Distance for each of the two stars.
    //    nodeOffset           siblingOffset
    //  v-------------v v----------------------v
    //  O ------------ x --------------------- o
    //  ^              ^                       ^
    // node (0)       com          sibling (maximumSeparation)

    float com = (sibling.nodeMass * separation) / (node.nodeMass + sibling.nodeMass);
    float nodeOffset = com;
    float siblingOffset = separation - com;

    // ...

Now, with the semimajor axes of both the node and its sibling determined, we can get to actually generating the orbital elements. However, there is one more thing we have to check: is this orbital node a star system or is it a barycenter?

In the abstract, this doesn’t matter; however, in our data structure for OrbitalNode, we have two pointers: starSystem and orbitalSystem. If this node represents a star system, then starSystem points to that object; if it represents a barycenter, then starSystem is null. When determining the orbit for a starSystem, we need to not only generate that star system’s orbitSystem, but reference that orbitSystem also in OrbitalNode.

    // If the node is a star system, update the star system's orbital system.
    if (node.starSystem != null)
        // If this is an L1 system, make sure the inclination is zero.
        if (node.depth == 1)
        {node.starSystem.orbitalSystem.GenerateRandom(nodeOffset, 0, 0.15f);}

        // Otherwise, allow inclination.
        {node.starSystem.orbitalSystem.GenerateRandom(nodeOffset, 0.15f);}

        // Set the orbital system.
        node.orbitalSystem = node.starSystem.orbitalSystem;
    // Otherwise this is a barycenter.
        node.orbitalSystem = new OrbitalSystem();

        // If this is an L1 system, make sure the inclination is zero.
        if (node.depth == 1)
        {node.orbitalSystem.GenerateRandom(nodeOffset, 0, 0.15f);}

        // Otherwise, allow inclination.
        {node.orbitalSystem.GenerateRandom(nodeOffset, 0.15f);}

    // Store the barycenter offset.
    node.barycenterSeparation = nodeOffset;

Generating the Elements

At last, we have finally come to OrbitalSystem.GenerateRandom(), where we actually hack out the orbital elements.

We have two overloads of GenerateRandom(); both take the semimajor axis and the deviation as an argument, but one also takes in a specified inclination as an argument.

The “deviation” argument is a measure of how exotic we want to allow our orbit to get, from 0 (circular with no inclination) to 1 (parabolic), where 1 is an exclusive boundary.

The reason we want to specify inclination sometimes is that for L1 pairs, we always set inclination to 0. This is because the inclination of the solar system at large, with respect to the galaxy, is accounted for in the SolarSystem’s OrbitalSystem, and need not be accounted for down in the StarSystem’s OrbitalSystems, which are all defined in the context of the SolarSystem’s plane.

With this, we can generate random orbital elements like so:

public void GenerateRandom(float smaj, float inc, float deviation) 
    // Determine Orbital Elements
    semimajorAxis = smaj;
    eccentricity = Numericals.randomHalfGaussian(0, deviation);
    inclination = inc;
    longitudeOfAscendingNode = Numericals.randomInt(0, 360);
    argumentOfPeriapsis = Numericals.randomInt(0, 360);

    // Determine other orbital information 
    semiminorAxis = CalculateSemiminor(semimajorAxis, eccentricity);
    semiLatusRectum = semimajorAxis * (float) (1 - Math.Pow(eccentricity, 2));

public void GenerateRandom(float smaj, float deviation)
    float randomInc = Numericals.randomHalfGaussian(0, (deviation * 45));
    GenerateRandom(smaj, randomInc, deviation);

Here we encounter another Numericals method: RandomHalfGaussian. Random Half Gaussian is exactly what it sounds like: it returns a normally distributed random value centered on the first argument (the mean), using the second argument (deviation), but it will only return numbers above the mean, never below it. Basically, imagine we have lopped off the left side of the gaussian distribution.

Determination of the Perifocal Coordinate System

It turns out that simply generating the 5 orbital elements is not sufficient information when it comes time to propagate these bodies or plot their elliptical points. In order to do this, we will also need the Perifocal Vectors of the orbital system.

The Perifocal Coordinate Frame of an orbit is a set of three orthogonal vectors, P, Q and W, which define a reference frame for an orbit, in which P points to the orbit’s periapsis, Q is 90 degrees past periapsis, and W is orthogonal to the orbital plane.

The perifocal reference frame will turn out to be an indispensable tool for us when it comes time to plot the elliptical points of orbits and to propagate orbital bodies.

Determining the reference frame is fairly simple, but an explanation of where these equations comes from exceeds the scope of this blog post. Suffice to say,

// Determine the perifocal vectors P Q and W for this orbit, as related to the equatorial coordinate system.
private void DeterminePerifocalVectorsFromOES()
    // Convert orbital elements of degrees into radians.
    float lanRadians = (float) ((Math.PI / 180) * longitudeOfAscendingNode);
    float aopRadians = (float) ((Math.PI / 180) * argumentOfPeriapsis);
    float incRadians = (float) ((Math.PI / 180) * inclination);

    // Determine P components
    float pi = (float) (Math.Cos(lanRadians) * Math.Cos(aopRadians)
        - Math.Sin(lanRadians) * Math.Cos(incRadians) * Math.Sin(aopRadians));
    float pj = (float) (Math.Sin(lanRadians) * Math.Cos(aopRadians)
        + Math.Cos(lanRadians) * Math.Cos(incRadians) * Math.Sin(aopRadians));
    float pk = (float) (Math.Sin(incRadians) * Math.Sin(aopRadians));

    // Determine Q components
    float qi = (float) (-Math.Cos(lanRadians) * Math.Sin(aopRadians)
        - Math.Sin(lanRadians) * Math.Cos(incRadians) * Math.Cos(aopRadians));
    float qj = (float) (-Math.Sin(lanRadians) * Math.Sin(aopRadians)
        + Math.Cos(lanRadians) * Math.Cos(incRadians) * Math.Cos(aopRadians));
    float qk = (float) (Math.Sin(incRadians) * Math.Cos(aopRadians));

    // Determine W components
    float wi = (float) (Math.Sin(incRadians) * Math.Sin(lanRadians));
    float wj = (float) (-Math.Sin(incRadians) * Math.Cos(lanRadians));
    float wk = (float) (Math.Cos(incRadians));

    // Assemble Vectors
    P = new System.Numerics.Vector3(pi, pj, pk);
    Q = new System.Numerics.Vector3(qi, qj, qk);
    W = new System.Numerics.Vector3(wi, wj, wk);

If you are interested in where these expressions came from, I refer you once again to Bate, Mueller and White’s textbook, “Fundamentals of Astrodynamics,” where I learned these expressions.

Generating Opposing Binary Orbits

So far, in our quest to generate orbital elements for a sibling pair, we have only generated the elements for node and not sibling. The Sibling’s orbit will be very similar to Node’s, save for two key differences:

  1. Sibling has a different semimajor axis
  2. Sibling’s Argument of Perigee will be offset by Pi (180 degrees)

And that’s actually it.

In GenerateSiblingKeplerianOrbits(),

    if (sibling.starSystem != null)
        sibling.starSystem.orbitalSystem.GenerateOpposingBinaryOrbit(siblingOffset, node.orbitalSystem);
        sibling.orbitalSystem = sibling.starSystem.orbitalSystem;
        sibling.orbitalSystem = new OrbitalSystem();
        sibling.orbitalSystem.GenerateOpposingBinaryOrbit(siblingOffset, node.orbitalSystem);

    // Store the barycenter offset.
    sibling.barycenterSeparation = siblingOffset;

Note that we are checking whether this is a barycenter or star system, just like earlier, with the preliminary node.

GenerateOpposingBinaryOrbit simply transcribes the relevant orbital elements while taking in the new semimajor axis. It also flips the argument of perigee by 180 degrees.

        public void GenerateOpposingBinaryOrbit(float smaj, OrbitalSystem opposingSystem)
            semimajorAxis = smaj;
            eccentricity = opposingSystem.eccentricity;
            inclination = opposingSystem.inclination;
            longitudeOfAscendingNode = opposingSystem.longitudeOfAscendingNode; 
            argumentOfPeriapsis = (opposingSystem.argumentOfPeriapsis + 180 > 360) ? opposingSystem.argumentOfPeriapsis - 180 : opposingSystem.argumentOfPeriapsis + 180;

            // Determine other orbital information 
            semiminorAxis = CalculateSemiminor(semimajorAxis, eccentricity);
            semiLatusRectum = semimajorAxis * (float) (1 - Math.Pow(eccentricity, 2));

Last, we need GenerateSiblingKeplerianOrbits() to recursively call itself.

        // Generate the Keplerian orbits for this node's children and its siblings children, if they have any.
        // Only call this on one child; the method is responsible for finding the child's sibling.
        if (node.child1 != null)
            generateSiblingKeplerianOrbits(node.child1, maximumDepth, maximumSeparation - (maximumSeparation / (maximumDepth+1)));

        if (sibling.child1 != null)
            generateSiblingKeplerianOrbits(sibling.child1, maximumDepth, maximumSeparation - (maximumSeparation / (maximumDepth+1)));

The Results

We have not yet explained how to visualize these results, but here are some in-software screenshots of the OrbitalSystems we just created.

Here we can see a binary solar system in which each star has a roughly (but not perfectly) circular orbit. The blue star (a type A star) is of greater mass than the yellow (type G) star. Between the two star systems, there are 23 planetary systems.

Roughly circular binary system, viewed from above

Here we can see a binary solar system in which the stars are in a much more elliptical set of orbits. Here it becomes evident that the purple orbit’s right focus is centered on the shared barycenter, while the yellow orbit’s left focus is centered on the shared barycenter.

A more eccentric binary system, viewed from above

Things get quite spicy when we add a third star system. In the next screenshot, you can see that the two suns on the right act like a single body, and that they, as a pair, are in a binary relationship with the star on the left.

You can see also that they, as a pair, have a much higher mass, and therefore are much closer to the center of the system.

A trinary system, viewed at a slight angle

Let’s turn things up with a quarternary system:

At this level, things get somewhat difficult to visually interpret. Let’s dive into what we’re seeing here:

In the next screen shot, you can see that the left binary pair has been encircled in blue. You can see that it is moving, as a pair, along the ellipse indicated by the blue vector. This ellipse is centered on the white X at the center of the shot.

The right binary pair is encircled in red, and moves along the ellipse indicated by the red vector. The ellipse is also centered on the white X at the center of the shot.

Next Steps

Next up, we’ll be visualizing these orbits in Unity. Like most things, visualizing the content turns out to be just as hard as generating the content.

After that, we’ll be propagating the bodies: that is, moving them along their orbits for the user to see.

After that, we’ll be making this Solar view interactable, so that when a user clicks on a star, we can jump to that particular star system and interact with its planets (which will be icospherical world models)!

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s