Help with Creating Orbital Manoeuvres

Unity Engine
, , , , ,
1

I am making a 3d orbital fleet combat game and trying to add a manoeuvre node system like seen in ksp. I have created a system that can define orbits and run the simulation as well as implement a timewarp system however turning deltaV into new orbits has proven dificult, especially with defining the Argument of Periapsis, Longitude of Ascending Node and the position vector. the for the orbital mechanics simulation is below:

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

[ExecuteInEditMode]
public class OrbitalMechanicsSimulation : MonoBehaviour
{
    public float altitudeKm;
    [Range(0.1f, float.PositiveInfinity)] public float semiMajorAxis;
    [Range(0, 0.9999999701976775f)] public float eccentricity;
    public float inclination; //Inclination in degrees
    public float argumentOfPeriapsis; //Argument of periapsis in degrees
    public float longitudeOfAscendingNode; //Longitude of ascending node in degrees
    public float orbitalPeriod; //Orbital period in seconds

    public Transform centrePoint;
    public float centralBodyMass = 1f; 
    public static float gravitationalConstant = 6.67430e-11f;

    #region Previous Values
    // Track previous values to detect changes
    private float previousSemiMajorAxis;
    private float previousEccentricity;
    private float previousInclination;
    private float previousArgumentOfPeriapsis;
    private float previousLongitudeOfAscendingNode;
    private float previousCentralBodyMass;
    private float previousTimewarpFactor;
    private float previousAltitudeKm;
    #endregion

    [SerializeField] private float currentTime = 0f;
    
    private LineRenderer lr;
    public int numPoints = 100;

    public float timewarpFactor = 1;
    [SerializeField] float quietTimeScaleFactor = 77284.3f;

    universalConstants _UniversalConstants;

    public Vector3 OrbitalManoeuvre;
    public Vector3 velocity;
    public Vector3 ManoeuvreNode;
    public float acceleration;
    public float burntime;

    public float orbitalSpeed;
    float trueAnomaly;
    [SerializeField] GameObject lrHolder;


    float deltaSemiMajorAxis;
    float deltaEccentricity;
    float deltaInclination;
    float deltaArgumentOfPeriapsis;
    float deltaLongitudeOfAscendingNode;
    float deltaTrueAnomaly;

    void Start()
    {
        if(centrePoint == null)
        {
            Debug.LogError("centrePoint is not set");
        }

        lr = GetComponent<LineRenderer>();
        if (lr != null)
        {
            lr.positionCount = numPoints + 1;
        }

        UpdatePreviousValues();
        UpdateOrbitalPeriod();

        _UniversalConstants = GameObject.Find("/GameMaster").GetComponent<universalConstants>();
        quietTimeScaleFactor = _UniversalConstants.QuietTimeScaleFactorStorage;
        centralBodyMass = centrePoint.GetComponent<planetaryStats>().BodyMass;
    }

    void Update()
    {
        orbitalSpeed= velocity.magnitude;

        timewarpFactor = _UniversalConstants.TimeWarpValue;
        currentTime += Time.deltaTime;
        if (previousAltitudeKm != altitudeKm)
        {
            semiMajorAxis = (centrePoint.localScale.x / 2) + (altitudeKm / 100);
            previousAltitudeKm = altitudeKm;    
        }

        if (NewOrbitalParameters())
        {
            UpdateOrbitalPeriod();
            UpdatePreviousValues();
        }


        if (orbitalPeriod <= 0f) { return; }
        
        if (Application.isPlaying)
        {
            RunSimulation();
        }
        
        DrawOrbit();

        if (Input.GetKeyDown(KeyCode.Space))
        {
            OrbitalManoeuvring();
        }
    }

    bool NewOrbitalParameters() { return semiMajorAxis != previousSemiMajorAxis || eccentricity != previousEccentricity || inclination != previousInclination || argumentOfPeriapsis != previousArgumentOfPeriapsis || longitudeOfAscendingNode != previousLongitudeOfAscendingNode || centralBodyMass != previousCentralBodyMass || timewarpFactor != previousTimewarpFactor; } //Check if any parameter has changed

    #region Updating orbital period
    void UpdateOrbitalPeriod()
    {
        float orbitalPercentage = 0f;
        if (orbitalPeriod != 0) { orbitalPercentage = currentTime / orbitalPeriod; }
        //Debug.Log(orbitalPercentage);
        orbitalPeriod = (2 * Mathf.PI * Mathf.Sqrt(Mathf.Pow(semiMajorAxis, 3) / (gravitationalConstant * centralBodyMass))) / (timewarpFactor * quietTimeScaleFactor);
        //Debug.Log(orbitalPeriod);
        if (orbitalPercentage * orbitalPeriod < float.MaxValue) { currentTime = orbitalPercentage * orbitalPeriod; }
        //Debug.Log(currentTime);
    }

    void UpdatePreviousValues()
    {
        previousSemiMajorAxis = semiMajorAxis;
        previousEccentricity = eccentricity;
        previousInclination = inclination;
        previousArgumentOfPeriapsis = argumentOfPeriapsis;
        previousLongitudeOfAscendingNode = longitudeOfAscendingNode;
        previousCentralBodyMass = centralBodyMass;
        previousTimewarpFactor = timewarpFactor;
    }
    #endregion

    void RunSimulation()
    {
        float meanAnomaly = 2 * Mathf.PI * currentTime / orbitalPeriod;
        float eccentricAnomaly = MeanToEccentricAnomaly(meanAnomaly);
        trueAnomaly = EccentricToTrueAnomaly(eccentricAnomaly);
        float r = semiMajorAxis * (1 - eccentricity * eccentricity) / (1 + eccentricity * Mathf.Cos(trueAnomaly));

        float orbitalSpeed = Mathf.Sqrt(gravitationalConstant * centralBodyMass / r);

        Vector3 orbitalPosition = new Vector3(r * Mathf.Cos(trueAnomaly), 0, r * Mathf.Sin(trueAnomaly));
        orbitalPosition = ApplyOrbitalTransform(orbitalPosition);

        transform.position = orbitalPosition + centrePoint.position;

        Vector3 direction = (transform.position - centrePoint.position).normalized;
        velocity = direction * orbitalSpeed * timewarpFactor * quietTimeScaleFactor;
    }

    #region RunSimulation extras
    float MeanToEccentricAnomaly(float meanAnomaly)
    {
        float eccentricAnomaly = meanAnomaly;
        for (int i = 0; i < 10; i++)
        {
            eccentricAnomaly = eccentricAnomaly - (eccentricAnomaly - eccentricity * Mathf.Sin(eccentricAnomaly) - meanAnomaly) / (1 - eccentricity * Mathf.Cos(eccentricAnomaly));
        }
        return eccentricAnomaly;
    }

    float EccentricToTrueAnomaly(float eccentricAnomaly)
    {
        return 2 * Mathf.Atan2(Mathf.Sqrt(1 + eccentricity) * Mathf.Sin(eccentricAnomaly / 2), Mathf.Sqrt(1 - eccentricity) * Mathf.Cos(eccentricAnomaly / 2));
    }

    Vector3 ApplyOrbitalTransform(Vector3 position)
    {
        Quaternion inclinationRotation = Quaternion.Euler(inclination, 0f, 0f);
        position = inclinationRotation * position;

        Quaternion periapsisRotation = Quaternion.Euler(0f, argumentOfPeriapsis, 0f);
        position = periapsisRotation * position;
        
        Quaternion ascendingNodeRotation = Quaternion.Euler(0f, longitudeOfAscendingNode, 0f);
        position = ascendingNodeRotation * position;

        return position;
    }
    #endregion
    
    void DrawOrbit()
    {
        if (lr == null) { return; }

        Vector3[] positions = new Vector3[numPoints + 1];
        float deltaTheta = 2 * Mathf.PI / numPoints;

        for (int i = 0; i < numPoints; i++)
        {
            float theta = i * deltaTheta;
            float r = semiMajorAxis * (1 - eccentricity * eccentricity) / (1 + eccentricity * Mathf.Cos(theta));
            float x = r * Mathf.Cos(theta);
            float z = r * Mathf.Sin(theta);

            Vector3 orbitalPosition = new Vector3(x, 0, z);

            orbitalPosition = ApplyOrbitalTransform(orbitalPosition); 

            orbitalPosition += centrePoint.position;
            
            if (i == 0)
            {
                positions[numPoints] = orbitalPosition;
            }
            positions[i] = orbitalPosition;
        }

        lr.SetPositions(positions);
    }
    


    void OrbitalManoeuvring()
    {
        Vector3 newVelocity = new Vector3(velocity.x * 100f, velocity.y * 100f, velocity.z * 100f);

        float gameSpeedUnitsToMetresPerSeconds = 10000;

        float deltaTangential = OrbitalManoeuvre.x / gameSpeedUnitsToMetresPerSeconds;
        float deltaOutOfPlane = OrbitalManoeuvre.y / gameSpeedUnitsToMetresPerSeconds;
        float deltaRadial = OrbitalManoeuvre.z / gameSpeedUnitsToMetresPerSeconds;

        Vector3 radialDirection = (transform.position - centrePoint.position).normalized;

        Vector3 orbitalNormal = Vector3.Cross(radialDirection, newVelocity).normalized;
        if (orbitalNormal == Vector3.zero)
        {
            if (radialDirection != Vector3.up)
                orbitalNormal = Vector3.Cross(radialDirection, Vector3.up).normalized;
            else
                orbitalNormal = Vector3.Cross(radialDirection, Vector3.right).normalized;
        }

        Vector3 tangentialDirection = newVelocity.normalized;

        float radial = Vector3.Dot(newVelocity, radialDirection);
        float tangential = Vector3.Dot(newVelocity, tangentialDirection);
        float outOfPlane = Vector3.Dot(newVelocity, orbitalNormal);

        float GM = gravitationalConstant * (centralBodyMass * 5.97e24f);
        float radialDistance = Vector3.Distance(this.transform.position, centrePoint.position);
        float specificAngularMomentum = radialDistance * velocity.magnitude;
        float perifocalDistance = semiMajorAxis * (1 - eccentricity * eccentricity);

        //Debug.Log("Perifocal Distance" + perifocalDistance);
        //Debug.Log("True Anomaly " + trueAnomaly);
        //Debug.Log("SpecificAngularMomentum" + specificAngularMomentum);

        //Debug.Log("Radial Direction: " + radialDirection);
        //Debug.Log("Orbital Normal: " + orbitalNormal);
        //Debug.Log("Tangential Direction: " + tangentialDirection);
        //Debug.Log("Velocity: " + newVelocity);

        //Debug.Log("radial: " + radial);
        //Debug.Log("tangential: " + tangential);
        //Debug.Log("outOfPlane: " + outOfPlane);

        deltaSemiMajorAxis = 2 * ((semiMajorAxis * semiMajorAxis) / GM) * ((radial / velocity.magnitude) * deltaRadial + (tangential / velocity.magnitude) * deltaTangential);
        deltaEccentricity = (1 / GM) * ((2 * (GM / radialDistance) - (GM / (2 * semiMajorAxis))) * (deltaTangential / tangential) - (1 - (radialDistance / semiMajorAxis)) * (tangential / GM) * deltaRadial);
        deltaInclination = ((radialDistance * Mathf.Cos(trueAnomaly + argumentOfPeriapsis)) / specificAngularMomentum) * deltaOutOfPlane;
        deltaArgumentOfPeriapsis = (((radialDistance * Mathf.Sin(trueAnomaly)) / (specificAngularMomentum * eccentricity)) * deltaRadial) - (((radialDistance * Mathf.Cos(trueAnomaly)) / (specificAngularMomentum * eccentricity)) * (1 + radialDistance / perifocalDistance) * deltaTangential) - ((Mathf.Sin(inclination)) / (specificAngularMomentum * Mathf.Sin(inclination))) * deltaOutOfPlane;
        deltaLongitudeOfAscendingNode = ((radialDistance * Mathf.Sin(trueAnomaly + argumentOfPeriapsis)) / (specificAngularMomentum * Mathf.Sin(inclination))) * deltaOutOfPlane;

        deltaTrueAnomaly = (specificAngularMomentum / (radialDistance * eccentricity)) * (-((radialDistance / perifocalDistance) * Mathf.Sin(trueAnomaly) * deltaRadial) + (1 + (radialDistance / perifocalDistance)) * Mathf.Cos(trueAnomaly) * deltaTangential);


        //        Debug.Log("delta Semi-Major Axis: " + deltaSemiMajorAxis);
        //        Debug.Log("delta Eccentricity: " + deltaEccentricity);
        //        Debug.Log("delta Inclination: " + deltaInclination);
        Debug.Log("delta Argument of Periapsis: " + deltaArgumentOfPeriapsis);
        Debug.Log("delta Longitude of Ascending Node: " + deltaLongitudeOfAscendingNode);
        //        Debug.Log("delta true anomaly: " + deltaTrueAnomaly);

        displayNewOrbit();
    }

    void displayNewOrbit()
    {
        GameObject newLrHolder = Instantiate(lrHolder);
        LineRenderer newLr = newLrHolder.GetComponent<LineRenderer>();

        if (newLr == null) { return; }
        newLr.positionCount = numPoints + 1;


        float newSemiMajorAxis = semiMajorAxis + deltaSemiMajorAxis;
        float newEccentricity = eccentricity + deltaEccentricity;
        //Debug.Log("New Semi-Major Axis: " + newSemiMajorAxis);
        //Debug.Log("New Eccentricity: " + newEccentricity);


        Vector3[] positions = new Vector3[numPoints + 1];
        float deltaTheta = 2 * Mathf.PI / numPoints;

        for (int i = 0; i < numPoints; i++)
        {
            float theta = i * deltaTheta;
            float r = newSemiMajorAxis * (1 - newEccentricity * newEccentricity) / (1 + newEccentricity * Mathf.Cos(theta));
            float x = r * Mathf.Cos(theta);
            float z = r * Mathf.Sin(theta);

            Vector3 orbitalPosition = new Vector3(x, 0, z);

            orbitalPosition = newApplyOrbitalTransform(orbitalPosition);

            orbitalPosition += centrePoint.position;

            if (i == 0)
            {
                positions[numPoints] = orbitalPosition;
            }
            positions[i] = orbitalPosition;
        }

        newLr.SetPositions(positions);
    }

    Vector3 newApplyOrbitalTransform(Vector3 position)
    {
        float newInclination = inclination + deltaInclination;
        float newArgumentOfPeriapsis = argumentOfPeriapsis + deltaArgumentOfPeriapsis;
        float newLongitudeOfAscendingNode = longitudeOfAscendingNode + deltaLongitudeOfAscendingNode;
        // Debug.Log("New Inclination: " + newInclination);
        //Debug.Log("New Argument of Periapsis: " + newArgumentOfPeriapsis);
        //Debug.Log("New Longitude of Ascending Node: " + newLongitudeOfAscendingNode);

        Quaternion inclinationRotation = Quaternion.Euler(newInclination, 0f, 0f);
        position = inclinationRotation * position;

        Quaternion periapsisRotation = Quaternion.Euler(0f, newArgumentOfPeriapsis, 0f);
        position = periapsisRotation * position;

        Quaternion ascendingNodeRotation = Quaternion.Euler(0f, newLongitudeOfAscendingNode, 0f);
        position = ascendingNodeRotation * position;

        return position;
    }

    void applyNewOrbit()
    {

    }

}

Any help would be greatly appreciated as I have completely stalled.

2

Woah, slow down … :smile:

That’s very technical. Can’t offer any direct help BUT to say that these sort of “simulations” often deviate dramatically from both the real-world and simulations intended to accurately simulate the real world.

KSP may be an exception, I recall they even made their own double-precision physics simulation. The kind of thing you will have to have for certain if the orbit around the planet exceeds 5 to 10 thousand units (assuming planet origin is at 0,0,0). Otherwise you’ll run into floating point accuracy issues, much like “Minecraft far travel” videos showcase.

It’s quite likely that you won’t get very far with those “universal constants” for the same reason as a float type may not be accurate enough to represent some of them and when you put them through even a simple calculation the result could already contain rounding errors, more so the bigger the number to the left of the decimal is. Even double types may have such issues considering the scale of something like a planet, let alone a solar system or even a galaxy.

The alternative is to not lean into the technical side but rather simplify. Games often employ tricks. Here I would tend to orbit in a circle at discrete distances - much like electrons orbit an atom at discrete energy levels. When changing orbit “levels” you can still smoothly animate that.

3

I’ll second CodeSmile’s “just make it fun” mantra.

I’ve not played KSP but I do know orbital mechanics and I also live in Earth’s gravity well. :slight_smile:

Looking at a video on maneuver nodes, if I had to guess, here’s several things going on:

  • a really amazing UI system that gives you fine-grain control of both vector and scalar direct and emergent values on a 2D screen
  • a solver to integrate results into the future and compute trajectories and display them
  • perhaps a custom high-precision math system, but likely double is adequate.

I can guarantee that float is inadequate to the task. Horribly inadequate in fact:

https://starmanta.gitbooks.io/unitytipsredux/content/floating-point.html

“Think of [floating point] as JPEG of numbers.” - orionsyndrome on the Unity3D Forums

Remember also that computer physics are always a discrete iterative simulation of what in reality is a continuous process. This always introduces an error term, and this error term almost always leads to either decay or instability of orbits, particularly when angular bearings between orbiting objects are changing rapidly per sample interval (eg, during close orbital passes).

At the end of the day you’re making a game, so keep in mind what might be fun. :slight_smile:

4

Not quite. :wink:

There’s the smallest unit of time, the Planck Time, which means our world does run on discrete time steps. However, we still have a problem with Heisenburg’s uncertainty principle plus all those curled up extra dimensions that may or may not exist.

The Planck time is about 5.391247 × 10^-44 seconds or veryveryveryveryveryveryveryveryveryvery short.

If I expanded it right, the number is:
0.000000000000000000000000000000000000000000005391247 seconds per universe tick

1 Like
5

For anyone wondering, I solved this by rewriting my code to define orbits based on Velocity and vector3 position then calculating with that.