Generating the Positions and Velocities of Orbital Bodies

This post builds on the concepts contained in “Generating Orbital Elements from Abstract Solar Hierarchies.”

Context

So far in our journey to generate and visualize the orbits of N-ary Simplex Solar Systems, we have:

  1. Generated the simplex solar hierarchy in the abstract, with no geometry.
  2. Translated this abstract hierarchy into a set of genuine Keplerian orbital elements.

Now, for our final step before visualization, we must actually generate a series of positions and velocities for these orbital bodies, so we can draw them to screen and animate them. At this point, we have all the data we need.


Data Structure

Our Positions and Velocities will each be represented by an array. Each index of the array corresponds with the same index in the opposite array. That is, for positions[n], the velocity at that position is represented by velocities[n].

We will have to decide on an orbital resolution – the number of entries in our position and velocity arrays. This value will of course be a configurable one. For our initial exploration, 360 seems like a nice resolution to choose. This means that positions and velocities will be arrays of size 360.

Each entry will contain a vector3; in the positions array, it will be a vector indicating the position of the planet relative to its orbital focus; in the velocities array, it will be a vector indicating the instantaneous velocity of the body at the corresponding position, relative to its perifocal coordinate frame. If you’re not sure what a perifocal coordinate frame is, we discuss it in a previous post.

Each entry in our position and velocity arrays represents the position-velocity at a specified true anomaly. What exactly true anomaly is, as opposed to mean anomaly and eccentric anomaly, is out of scope for this post. Suffice to say, true anomaly is the angle between the orbital body’s periapsis point, the orbital focus, and the orbital body. It is the rotational angle of the body as viewed from the orbital focus.


The Position Array

The Mathematics

To determine the position of the orbital body at a given true anomaly, we use the expression:

Where:

  • r is the position vector relative to the orbital focus.
  • r is the conic section.
  • v is the true anomaly.
  • P, Q are the orthogonal components of the parent perifocal coordinate frame.

The Conic Section

The conic section is the curve obtained by intersecting a cone and a plane. The curves obtained this way may be circles, ellipses, parabolas and hyperbolas. As it turns out, any orbit can be described by a conic section, if we use a cone of the proper angle and a plane with the correct orientation.

Look Familiar?

The equation for the conic section is as follows:

Where:

  • r is the distance from the orbital focus at the specified anomaly
  • p is the semi-latus rectum
  • v is the true anomaly.

The semi-latus rectum is simply the distance of a chord drawn from an orbital focus to the edge of the ellipse, orthogonal to the semimajor axis. You can google how to calculate it; the recursive equation train stops here.

The Code

With the expression for position, we are ready to create a method to calculate the position of a body given an anomaly.

// Determine the position vector of a body along the ellipse, given an anomaly and a specific set of perifocal vectors.
public Vector3 DeterminePointAtAnomaly(float anomalyDeg, Vector3 p, Vector3 q)
{
    // Convert anomaly from degrees to radians
    float anomalyRad = (float) ((Math.PI / 180) * anomalyDeg); 

    //  Calculate Conic
    float r = CalculateConicAtAnomaly(anomalyRad);

    // Return the relative position
    return (r * (float) Math.Cos(anomalyRad) * p) + (r * (float) Math.Sin(anomalyRad) * q); 
}

This method makes use of another method: CalculateConicAtAnomaly.

// Evaluate the Polar Equation of the Conic Section at a particular angle.
private float CalculateConicAtAnomaly(float anomalyRadians)
{return semiLatusRectum / (1 + eccentricity * (float) Math.Cos(anomalyRadians));}

Invoking DeterminePointAtAnomaly() for each point along the orbital resolution, we are able to populate the position array.


The Velocity Array

The Mathematics

To determine the velocity of the orbital body at a given true anomaly, we use the expression:

Where:

  • v = the velocity vector, relative to the orbital focus.
  • mu is the Standard Gravitational Parameter.
  • p is the semi-latus rectum.
  • P, Q are the orthogonal vectors of the perifocal reference frame.
  • v is the true anomaly.
  • e is the eccentricity.

The Standard Gravitational Parameter

The Standard Gravitational Parameter is a convenient shorthand for GM (The Gravitational Constant and the Body Mass). Mu is used to wrap these values together because the appear so often beside one another. In our solar system, the standard gravitational parameter of various celestial bodies is often known with more precision than the gravitational constant itself. More on that in another post.

In any case, for a binary star system with elliptical orbits, the SGP can be expressed:

Where:

  • a is the semimajor axis
  • T is the period of the orbit

The units of the SGP are often expressed in meters ^3, seconds ^-2. For us, in the timescales of orbiting star systems, AU/Years seems more convenient. We can of course convert our SGPs whenever we need to.

The Period

The period, of course, is the length of time that it takes for an orbital body to complete one full revolution around its orbital focus. The period can be expressed:

Where:

  • T is the period, in seconds
  • a-sum is the sum of the binary pair’s semimajor axes
  • G is the Gravitational Constant, 6.67430(15)×10−11 m3kg–1s–2
  • M1 and M2 are the masses of the two bodies

The Code

Equipped with these equations, we are now ready to determine the velocity at a given anomaly.

You will notice that in this method, I have separated out a number of different components of the velocity equation into separate variables. This is primarily because the velocity method was quite fucked for some time, and to debug it more easily, I pulled out each component of the equation so I could read its value.

When writing code for layered mathematical expressions I often find it cleaner and less error-prone to pull the equation apart into sub expressions. No programmer wants to see a 300 character-long line of hideous C# system method calls, it’s completely unreadable.

Without further ado:

// Determine the velocity vector of a body along the ellipse, given an anomaly and a specific set of perifocal vectors.
// This velocity is in AU per Year.
public void DetermineVelocityAtAnomaly(float anomaly, Vector3 p, Vector3 q)
{
    // Obtain Anomaly in Radians
    float anomalyRadians = (float) ((Math.PI / 180) * anomaly); 

    // Determine the square root of the standard gravitational parameter divided by the semi-latus rectum.
    double sqrtSgpOverSlr = Math.Sqrt(standardGravitationalParameter / semiLatusRectum);

    // Determine the P component of the velocity.
    Vector3 pComponent = -(Vector3.Multiply((float) Math.Sin(anomalyRadians), p));

    // Determine the Q component of the velocity.
    Vector3 qComponent = Vector3.Multiply((float) (eccentricity + Math.Cos(anomalyRadians)), q);

    // Sum the vector components.
    Vector3 pq = Vector3.Add(pComponent, qComponent);

    // Multiply by the square root of the standard gravitational parameter divided by the semi latus rectum.
    Vector3 velocity = Vector3.Multiply((float) sqrtSgpOverSlr, pq);

    // Update Collections.
    velocityVectors.Add(velocity);
    velocityMagnitudes.Add(velocity.Length());
}

Notice we are curating not only a velocity array, but a velocity magnitudes array. More on this later.

You will notice that the SGP is already evaluated before this method is called. That’s because the calculation of the SGP and the Period are now being done in the generateRandom() methods of the OrbitSystem, where this is all taking place. To refresh yourself on OrbitSystem see this previous post.

In any case, here are the methods to determine the period and the SGP, in turn:

// Calculate the period of this two body system, in seconds.
public static float CalculatePeriod(float thisSemimajorAxisAU, float opposingSemimajorAxisAU, float thisMassSM, float opposingMassSM) 
{
    // Convert SMA from AU to M
    double thisSemimajorAxisM = (double) thisSemimajorAxisAU * Conversions.AU_TO_KM * 1000;
    double opposingSemimajorAxisM = (double) opposingSemimajorAxisAU * Conversions.AU_TO_KM * 1000;

    // Convert masses from AU to KG
    double thisMassKG = thisMassSM * Conversions.SM_TO_KG;
    double opposingMassKG = opposingMassSM * Conversions.SM_TO_KG;

    // Determine the period, in seconds.
    return (float) (2 * Math.PI * Math.Sqrt(
        Math.Pow((thisSemimajorAxisM + opposingSemimajorAxisM), 3) / (Constants.G * ((double) thisMassKG + (double) opposingMassKG))
    ));
}
// Calculate the Standard Gravitational Parameter for this system, in AU^3 / Years^2.
public static float CalculateStandardGravitationalParameterAUYears(float thisSemimajorAxisAU, float opposingSemimajorAxisAU, float periodSeconds)
{
    // float summedSmaAU = thisSemimajorAxisAU + opposingSemimajorAxisAU;
    float periodYears = periodSeconds / Conversions.SEC_PER_YR;

    // return (float) ((4 * Math.Pow(Math.PI, 2) * Math.Pow(summedSmaAU, 3)) / (Math.Pow(periodYears, 2)));
    return (float) ((4 * Math.Pow(Math.PI, 2) * Math.Pow(thisSemimajorAxisAU, 3)) / (Math.Pow(periodYears, 2)));
}

With these three methods, we are able to calculate the velocity vector and magnitude of the orbital body at any anomaly.


The Results

Pulling it All Together

Here we have built up two top-level methods: DeterminePointAtAnomaly() and DetermineVelocityAtAnomaly(). These two methods will need to be called in a coordinated way in order to populate the position[], velocityVector[] and velocityMagnitude[] arrays, in accordance with the orbital resolution. This is a straightforward and trivial task, which I leave in the hands of the capable reader.

What We Expect to See

When applying these positions and velocities to our game objects in unity (there will be an entire post dedicated to this), we hope to see a few important behaviors which indicate, intuitively, that our positions and velocities are really working and are mostly accurate.

  1. All binary pairs should always be at opposing anomalies.
  2. Bodies should be moving more swiftly at periapsis than at apoapsis.
A satellite, moving swiftly at periapsis and slowly at apoapsis, sweeping out equal areas of its ellipse in equal times.

Software Footage

Here is a video demonstration showing what our position and velocity arrays look like when propagated.

In this video we can see three demonstrations:

  1. A binary pair with roughly circular orbits.
  2. A binary pair with highly elliptcal orbits.
  3. A trinary solar system.

Generally, everything looks good. In the highly elliptical system, we are able to perceive the difference in speed at perigee and apogee, as we had hoped.

The working Trinary system lets us know that our propagation holds for all layers of the orbital hierarchy, and that our relative, tiered perifocal frames are all communicating.


Known Issues; Next Steps

In the video above, two issues which become clear:

  1. In small-axis, high-speed orbits (small, highly elliptical orbits), star positions start out correct but begin to drift within a few revs (see the left pair in the Trinary system above). I believe I know the source of this problem, and it is not in the position or velocity calculations, but in the Propagator, which we haven’t discussed yet.
  2. Something is occasionally goofy with the planetary orbits. The circular speed of far-off planets should be notably less than that of near-planets; sometimes, it isn’t. Probably this has to do with some error numerical or rounding resulting from the fact that unlike stars, a planet’s mass is negligible compared to that of a star. Instead of using a two body orbital system for the planets, I probably will need to switch to a point-satellite orbital system. More on this later.

Next, we will figure out how to transform these position and velocity arrays into moving objects in Unity.

Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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