​

Wednesday, June 16, 2021

Gravity Part 2: Orbit Predictions (One-Body problem)

Introduction

This is the second post in the series I am creating, based on the logic from Physics-based Solar System Lessons:
The previous post sets the foundations in which part 2 and 3 are built upon. Part 1 was mandatory as there have been many changes since the original Physics-based Solar System Lessons post.

Gravity Part 2: Orbit Predictions (One-Body problem)

Once we have a system affected only by gravitational forces, the next logical step was to predict orbits of objects. This is not a trivial task, especially when there are multiple bodies in play. I focused on implementing object trajectory predictions for One-Body problems.

Orbital mechanics sketch (image credit: wikipedia [4])


Background

All the math and concepts behind the One/N-Body problems have been obtained from the books by Keeton, C. (2014) [1] and Walter, U. (2018) [2].

What is a One-Body problem? 
A system in which the source of gravity is stationary while another object is in motion (and therefore affected by the stationary object’s gravity).
 
What is a Two-Body Problem?
A system in which two objects interact with each other through their force of gravity.

Can we predict orbits?

There are 3 different problems to predicting trajectories:
  • One-Body problem: Solving the prediction is not complicated. The visual representation is simple as the trajectory will be fixed in time.
  • Two-Body problem: Things become tricky as both objects move in time. Any visualization would be complicated as the forces applied will change as objects move in time.
  • N-Body problem: This is a generalization of the Two-body problem to a larger scale. As mentioned before, visualization of the orbits would be complicated.
There are two possibilities to visualize the trajectories:
  • Given a system, let it run for a certain amount of time and we can visualize the position and velocity of the objects. This can be used to predict N-Body problems and see the future trajectories, but it will not give you any past trajectories.
  • Given a function obtain the past and future trajectory of an object. This is simple and doable for a One-Body problem.
Disclaimer: There are ways to solve N-Body problems numerically, but I have not tried to apply the math in Unity.

The first option can be used by setting up a system with the information provided in Gravity Part 1: An update. The second option is what I am going to explain in this post.

Predictions

As we are focused on One-Body problems, we will always have one of these trajectories:
  • Elliptical/Circular: object “captured” by the gravity of the massive object.
  • Parabolic: Object that is passing just fast enough to escape gravitation pull.
  • Hyperbolic: Object that is passing too fast and escapes gravitation pull.
Calculations for elliptical trajectories are easy, but those for parabolic and hyperbolic are a bit more complex. All these trajectories come from conic sections. The parabolic trajectory is a fringe case that I will not implement, but it can be added following the same logic as it is for elliptic/hyperbolic.

To create these conics I will be using the following data:
  • A Massive object that is the primary force attraction:
    • Given its position.
    • Given its gravity force.
  • Another object that we want to check its trajectory:
    • Given its position.
    • Given its velocity and direction.

The Code

From this point onwards I will be explaining the math as well as the code I have used. The code can be summarized into 2 parts:
  • Conics: Conic shape trajectories.
  • Orbit: Calculation of the trajectory using conic shapes.

Conics

As can be observed by the following picture the conic family can create the desired trajectories: Circular, Elliptical, parabolic, and hyperbolic.

Conic sections and orbits (image credit: braeunig [3])

All conics have the same basic structure, as such I created an interface IConic that affects all.
    
    
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

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

interface IConicSection
{
    double SemiMajorAxis { get; }
    double SemiMinorAxis { get; }
    double Eccentricity { get; }
    double Distance { get; }
    ReadOnlyCollection<Vector3> Points { get; }

    List<Vector3> GetPointsInWorldCoordinates(Transform objectTransform);
    void SetConicSectionTransform(Transform orbitPlane, Vector3 focusPoint, Vector3 pointInTrajectory, bool isVelocityIncreasing);
}


It is worth mentioning that both ellipse and hyperbola will get all equation variables in the constructor from the eccentricity and semi major axis. Here is an explanation of the public functions:

GetPointsInWorldCoordinates(...): 
Given a transform (rotation and traslation) obtain the trajectory points in world coordinates

SetConicSectionTransform(...):
Modify the Orbit Object transform component to be aligned with the conic section. We need to know whether the velocity is increasing to obtain the correct angle of the object in the conic segment.

Ellipse trajectory

This also considers the Circular trajectory as it is a fringe case of the ellipse. We need to solve the basic (no translation nor rotation) ellipse equation:

\[ \frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} = 1\]
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
using System;
using System.Collections.Generic;
using System.Collections.ObjectModel;
using System.Linq;
using UnityEngine;

public class Ellipse : IConicSection
{
    #region Properties
    public double SemiMajorAxis { get; private set; }
    public double SemiMinorAxis { get; private set; }
    public double Eccentricity { get; private set; }
    public double Distance { get; private set; } //Foci Distance
    public ReadOnlyCollection<Vector3> Points { get; private set; }
    #endregion

    #region Public Methods
    /// <summary>
    /// Basic ellipse in Unity default units, centered in the origin and no rotation applied
    ///   x^2/a^2 + y^2/b^2 = 1
    /// 
    ///   a = SemiMajorAxis
    ///   b = SemiMinorAxis
    /// </summary>
    /// <param name="eccentricity"> Ellipse eccentricity [0,1] </param>
    /// <param name="semiMajorAxis"> Length of the semi major axis </param>
    public Ellipse(double eccentricity, double semiMajorAxis, int numberOfPoints = 10)
    {
        Eccentricity = eccentricity;
        SemiMajorAxis = semiMajorAxis;
        SemiMinorAxis = CalculateSemiMinorAxis();
        Distance = CalculateFociDistance();

        CreateEllipseLocalPoints(numberOfPoints);
    }

    /// <summary>
    /// Get all the local points in world coordinates, applying:
    ///  - A specific transform (position, normal and rotation around local Y)
    /// </summary>
    /// <param name="objectTransform"> An object's transform </param>
    /// <returns> Ellipse points in world coordinates </returns>
    public List<Vector3> GetPointsInWorldCoordinates(Transform objectTransform)
    {
        return Points.Select(p => objectTransform.TransformPoint(p)).ToList();
    }

    /// <summary>
    /// Given the orbit plane, an ellipse focus point, an object's position in the 
    /// ellipse (and if it is accelerating), adjust the orbit plane position and its 
    /// rotation so a projection of the ellipse is correct in world coordinates.
    /// </summary>
    /// <param name="orbitPlane"> where the Ellipse will be projected </param>
    /// <param name="focusPoint"> ellipse focus </param>
    /// <param name="pointInTrajectory"> point in the ellipse edge </param>
    /// <param name="isVelocityIncreasing"> is point in trajectory increasing speed? used to get the correct angle </param>
    public void SetConicSectionTransform(Transform orbitPlane, Vector3 focusPoint, Vector3 pointInTrajectory, bool isVelocityIncreasing)
    {
        //Set position
        var distanceFromObject = Vector3.Distance(pointInTrajectory, focusPoint);
        var angleRadians = GetAngleBetweenPointInEllipseAndSemiMajorAxis(distanceFromObject, isVelocityIncreasing);

        var rot = Quaternion.AngleAxis(-(float)angleRadians * Mathf.Rad2Deg, orbitPlane.up);
        var semiMajorAxisDirection = (rot * (focusPoint - pointInTrajectory)).normalized;
        orbitPlane.position = focusPoint + semiMajorAxisDirection * (float)Distance;

        //Set rotation
        var semiMajorAxisDirectionLocal = orbitPlane.InverseTransformDirection(semiMajorAxisDirection);
        var angle = Vector3.Angle(semiMajorAxisDirectionLocal, Vector3.right);

        orbitPlane.Rotate(Vector3.up, angle, Space.Self);
        var localPos = orbitPlane.InverseTransformPoint(pointInTrajectory);

        if (!IsPointInEdge(localPos))
        {
            orbitPlane.Rotate(Vector3.up, -angle * 2, Space.Self);
        }
    }
    #endregion


    #region Private Methods
    private double CalculateSemiMinorAxis()
    {
        return SemiMajorAxis * Math.Sqrt(1 - Math.Pow(Eccentricity, 2));
    }

    private double CalculateFociDistance()
    {
        return Math.Sqrt(Math.Pow(SemiMajorAxis, 2) - Math.Pow(SemiMinorAxis, 2));
    }

    private void CreateEllipseLocalPoints(int numberOfPoints)
    {
        var points = new List<Vector3>();
        var gap = 360.0f / numberOfPoints;

        for (int i = 0; i < numberOfPoints; i++)
        {
            var theta = (i * gap) * Mathf.Deg2Rad;
            var x = SemiMajorAxis * Mathf.Cos(theta);
            var z = SemiMinorAxis * Mathf.Sin(theta);

            points.Add(new Vector3((float)x, 0, (float)z));
        }

        points.Add(points[0]);
        Points = points.AsReadOnly();
    }

    //Get the angle using the orbit equation
    private double GetAngleBetweenPointInEllipseAndSemiMajorAxis(double distance, bool isIncreasingVelocity)
    {
        var cosE = GetCenterCosAngleFromOrbitEquation(distance);
        var angle = Math.Acos((cosE - Eccentricity) / (1 - Eccentricity * cosE));

        return isIncreasingVelocity ? angle : -angle;
    }

    private double GetCenterCosAngleFromOrbitEquation(double distance)
    {
        return -(distance / SemiMajorAxis - 1) / Eccentricity;
    }

    /// <summary>
    /// Verify if the given point is in the ellipse
    /// </summary>
    /// <param name="point"> a point in 3D, note the Y-axis is not used here </param>
    /// <returns></returns>
    private bool IsPointInEdge(Vector3 point)
    {
        return Math.Abs(Math.Pow(point.x, 2) / Math.Pow(SemiMajorAxis, 2) + Math.Pow(point.z, 2) / Math.Pow(SemiMinorAxis, 2) - 1) < Constants.EPSILON;
    }
    #endregion
}

Note: to obtain the rotation with respect to the focus we use the orbit equation and its derivation. To follow the math I followed the formulas provided by Walter (equations 7:4:14) [2]
 

Hyperbolic trajectory:

We need to solve the basic (no translation nor rotation) hyperbolic equation:

\[ \frac{x^{2}}{a^{2}} - \frac{y^{2}}{b^{2}} = 1\]
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
using System;
using System.Collections.Generic;
using System.Collections.ObjectModel;
using System.Linq;
using UnityEngine;

public class Hyperbola : IConicSection
{
    #region Properties
    public double SemiMajorAxis { get; private set; }
    public double SemiMinorAxis { get; private set; }
    public double Eccentricity { get; private set; }
    public double Distance { get; private set; } //Minimum distance from foci to edge
    public double SemiLatusRectum { get; private set; }
    public ReadOnlyCollection<Vector3> Points { get; private set; }
    #endregion

    #region Public Methods
    /// <summary>
    /// Basic hyperbola in Unity default units, centered in the origin and no rotation applied
    ///   x^2/a^2 - y^2/b^2 = 1
    /// 
    ///   a = SemiMajorAxis
    ///   b = SemiMinorAxis
    /// </summary>
    /// <param name="eccentricity"> Hyperbola eccentricity > 1 </param>
    /// <param name="semiMajorAxis"> Length of the semi major axis </param>
    public Hyperbola(double eccentricity, double semiMajorAxis, int numberOfPoints = 10)
    {
        Eccentricity = eccentricity;
        SemiMajorAxis = semiMajorAxis;
        SemiMinorAxis = CalculateSemiMinorAxis();
        Distance = CalculateMinimumDistance();
        SemiLatusRectum = SemiMajorAxis * (1 - Math.Pow(Eccentricity, 2));        

        CreateHyperbolaLocalPoints(numberOfPoints);
    }

    /// <summary>
    /// Get all the local points in world coordinates, applying:
    ///  - A specific transform (position, normal and rotation around local Y)
    /// </summary>
    /// <param name="objectTransform"> An object's transform </param>
    /// <returns> Hyperbola points in world coordinates </returns>
    public List<Vector3> GetPointsInWorldCoordinates(Transform objectTransform)
    {
        return Points.Select(p => objectTransform.TransformPoint(p)).ToList();
    }

    /// <summary>
    /// Given the orbit plane, a hyperbola focus point, an object's position in the 
    /// hyperbola (and if it is accelerating), adjust the orbit plane position and its 
    /// rotation so a projection of the hyperbola is correct in world coordinates.
    /// </summary>
    /// <param name="orbitPlane"> where the Hyperbola will be projected </param>
    /// <param name="focusPoint"> hyperbola focus </param>
    /// <param name="pointInTrajectory"> point in the hyperbola edge </param>
    /// <param name="isVelocityIncreasing"> is point in trajectory increasing speed? </param>
    public void SetConicSectionTransform(Transform orbitPlane, Vector3 focusPoint, Vector3 pointInTrajectory, bool isVelocityIncreasing)
    {
        //Set position
        orbitPlane.position = focusPoint;

        //Set rotation
        float angle = GetAngle(orbitPlane, focusPoint, pointInTrajectory, isVelocityIncreasing);

        orbitPlane.Rotate(Vector3.up, angle, Space.Self);        
        var localPos = orbitPlane.InverseTransformPoint(pointInTrajectory);

        if (!IsPointInEdge(localPos))
        {
            orbitPlane.Rotate(Vector3.up, -angle * 2, Space.Self);
        }
    }
    #endregion


    #region Private Methods
    private double CalculateSemiMinorAxis()
    {
        return -SemiMajorAxis * Math.Sqrt(Math.Pow(Eccentricity, 2) - 1);
    }

    private double CalculateMinimumDistance()
    {
        return SemiMajorAxis * (1 - Eccentricity);
    }

    private void CreateHyperbolaLocalPoints(int numberOfPoints)
    {
        var pointsPositive = new List<Vector3>();
        var pointsNegative = new List<Vector3>();
        var offset = -SemiMajorAxis + Distance;

        for (int i = 0; i < numberOfPoints / 2; ++i)
        {
            var z = i * i / 5;
            var x = (float)(SemiMajorAxis * Math.Sqrt(1 + Math.Pow(z, 2) / Math.Pow(SemiMinorAxis, 2)) + offset);

            pointsPositive.Add(new Vector3(x, 0, z));

            if (i != 0)
            {
                pointsNegative.Add(new Vector3(x, 0, -z));
            }

        }

        pointsPositive.Reverse();
        pointsPositive.AddRange(pointsNegative);
        Points = pointsPositive.AsReadOnly();
    }

    private float GetAngle(Transform orbitPlane, Vector3 focusPoint, Vector3 pointInTrajectory, bool isVelocityIncreasing)
    {
        var distanceFromObject = Vector3.Distance(focusPoint, pointInTrajectory);
        var angleRadians = Math.Acos((SemiLatusRectum / distanceFromObject - 1) / Eccentricity);

        angleRadians = (isVelocityIncreasing) ? -angleRadians : angleRadians;

        var rot = Quaternion.AngleAxis((float)angleRadians * Mathf.Rad2Deg, orbitPlane.up);

        var semiMajorAxisDirection = (rot * (pointInTrajectory - focusPoint)).normalized;
        var semiMajorAxisDirectionLocal = orbitPlane.InverseTransformDirection(semiMajorAxisDirection);
        var angleLocal = Vector3.Angle(semiMajorAxisDirectionLocal, Vector3.right);

        return angleLocal;
    }

    /// <summary>
    /// Verify if the given point is in the hyperbola. Offset applied as the points are not
    /// centered in the "center" of the hyperbola, but 1 of the focus points.
    /// </summary>
    /// <param name="point"> a point in 3D, not the Y-axis is not used here </param>
    /// <returns></returns>
    private bool IsPointInEdge(Vector3 point)
    {
        var offset = new Vector3((float)(Distance - SemiMajorAxis), 0, 0);
        point -= offset;

        var distanceToEdge = Math.Abs(Math.Pow(point.x, 2) / Math.Pow(SemiMajorAxis, 2) - Math.Pow(point.z, 2) / Math.Pow(SemiMinorAxis, 2) - 1);
        return distanceToEdge < Constants.EPSILON;
    }
    #endregion
}

Note: The function GetAngle(...) returns the angle from Vector3.Angle which outputs a value [0, 180]. As it is missing the [180, 360] range we check if the object is in the trajectory.

Orbit

The Orbit code basically consists of the functionality that will trigger the orbit creation and visualization. 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
using System;
using UnityEngine;
using static Constants;

public class Orbit : MonoBehaviour
{
    public GameObject OrbitingObject { get; private set; }
    public GameObject OrbitAroundObject { get; private set; }
    public OrbitTrajectory OrbitTrajectory { get; private set; } = OrbitTrajectory.Undefined;

    private GameManager _gameManager;
    private LineRenderer _lineRenderer;
    private IConicSection _conicSection;

    public void CreateOrbit(GameObject orbitAroundObject, GameObject orbitingObject)
    {
        _gameManager = GameObject.Find("GameManager").GetComponent<GameManager>();

        OrbitingObject = orbitingObject;
        OrbitAroundObject = orbitAroundObject;

        var isVelocityIncreasing = OrbitingObject.GetComponent<SpaceObject>().IsVelocityIncreasing;

        _conicSection = CreateConicSection();
        _conicSection.SetConicSectionTransform(transform, OrbitAroundObject.transform.position, OrbitingObject.transform.position, isVelocityIncreasing);
        var realWorlTrajectoryPoints = _conicSection.GetPointsInWorldCoordinates(transform);

        _lineRenderer = GetComponent<LineRenderer>();
        _lineRenderer.positionCount = realWorlTrajectoryPoints.Count;
        _lineRenderer.SetPositions(realWorlTrajectoryPoints.ToArray());
    }

    private IConicSection CreateConicSection()
    {
        var focusPoint = OrbitAroundObject.transform.position;
        var point = OrbitingObject.transform.position;
        var distance = Vector3.Distance(point, focusPoint);

        var velocity = OrbitingObject.GetComponent<Rigidbody>().velocity;
        var mu = _gameManager.GravitationalConstantUnityScaled * OrbitAroundObject.GetComponent<SpaceObject>().Mass;

        var semiMajorAxis = Math.Pow(2 / (double)distance - Math.Pow(velocity.magnitude, 2) / mu, -1);
        var orbitNormal = -Vector3.Cross(velocity, focusPoint - point);

        transform.up = orbitNormal.normalized;

        var angularMomentum = orbitNormal.magnitude;
        var semiLatusRectum = Math.Pow(angularMomentum, 2) / mu;
        var eccentricity = Math.Sqrt(1 - semiLatusRectum / semiMajorAxis);
        
        IConicSection conicSection;

        if (IsObjectInEllipticalOrbit())
        {
            conicSection = new Ellipse(eccentricity, semiMajorAxis, 50);
            OrbitTrajectory = OrbitTrajectory.Elliptical;
        }
        else
        {
            conicSection = new Hyperbola(eccentricity, semiMajorAxis, 50);
            OrbitTrajectory = OrbitTrajectory.Escape;
        }

        return conicSection;
    }

    private bool IsObjectInEllipticalOrbit()
    {
        var ellipsePointVelocity = OrbitingObject.GetComponent<Rigidbody>().velocity;
        var distance = (OrbitAroundObject.transform.position - OrbitingObject.transform.position).magnitude;
        var escapeVelocity = OrbitAroundObject.GetComponent<SpaceObject>().GetEscapeVelocity(distance);

        return ellipsePointVelocity.magnitude < escapeVelocity;
    }
}

It is worth mentioning that all the conics need the eccentricity and semimajor axis to be calculated. These are obtained through simple formulas:

Semi major axis:
\[ a = \frac{1}{\sqrt{\frac{2}{d} -  \frac{v^{2}}{\mu} }}\]
\[ \mu  = G m\]
where d is the distance between objects, v is the velocity of the orbiting object, G is the gravitational constant and m is the mass of the gravity source.

Eccentricity:
\[ e  = \sqrt{1 - \frac{L}{a}}\]
\[ L  = \frac{p^{2}}{\mu}\]
where L is the Semi latus Rectum and p is the angular momentum.

One last thing worth mentioning is how I obtain the orbit normal vector. The code is a simplification of the actual operation. For better visibility I will put here the logic used:

1
2
3
4
Vector3 pointA = point + velocity;
Vector3 pointB = focusPoint;
Vector3 orbitNormal = Vector3.Cross(pointA - point, pointB - point);
//Vector3 orbitNormal = Vector3.Cross(-velocity, focusPoint - point);

Updates

After creating the orbit predictions, some changes were made to the Constants.cs and SpaceObject.cs files. 
 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
public static class Constants
{
    public static readonly double GRAVITATIONAL_CONSTANT_KM = 6.67430e-20; //N*km^2/kg^2
    public static readonly double EPSILON = 0.01;

    public const int FREE_VELOCITY = 0;
    public const int CIRCULAR_ORBIT_VELOCITY = 1;
    public const int ESCAPE_ORBIT_VELOCITY = 2;

    public const int UNDEFINED = -1;
    public const int COLLISION = 0;
    public const int ELLIPTICAL = 1;
    public const int ESCAPE = 2;    

    public enum InitialVelocity
    {
        Free = FREE_VELOCITY,
        CircularOrbit = CIRCULAR_ORBIT_VELOCITY,
        EscapeOrbit = ESCAPE_ORBIT_VELOCITY
    }

    public enum OrbitTrajectory
    {
        Undefined = UNDEFINED,
        Collision = COLLISION,
        Elliptical = ELLIPTICAL,
        Escape = ESCAPE
    }
}

In the constants I just created some variables to keep track of the state, though there is no functionality attached to it.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
using System;
using System.Linq;
using UnityEngine;
using static Constants;

public class SpaceObject : Gravity
{
    #region variables
    [Header("Orbit related")]
    [SerializeField] private GameObject _orbit;

    [Header("Velocity settings")]
    [SerializeField] private Vector3 _initialDirection = new Vector3(0, 0, 1);
    [SerializeField] private InitialVelocity _velocityType = InitialVelocity.Free;
    [SerializeField] private double _initialVelocity = 0; // In standard km/s
    
    [Header("SpaceObject Properties")]

    [SerializeField]
    [Tooltip("external gravity forces do not affect this object")]
    private bool _disableGravityAffected = false;

    [SerializeField]
    [Tooltip("Kilograms")]
    private double _mass = 0; //Kg

    [SerializeField]
    [Tooltip("Kilometers")]
    private double _radius = 0; //km

    [SerializeField]
    [Tooltip("Kilometers. Distance at which entry to atmosphere happens")]
    private double _entryInterfaceDistance = 0; //Km. distance at which entry to atmosphere happens

    private Rigidbody _rigidbody;
    private GameObject _orbitAround;
    private bool initializationComplete = false;

    [Header("Debug")]
    [SerializeField] private double _currentVelocity = 0; // In standard km/s
    #endregion

    #region Properties
    public double Mass => _mass;

    public double CollisionDistance { get => (_radius + _entryInterfaceDistance) / _gameManager.SpaceScaleKm; }

    public double CurrentVelocity
    {
        get { return _currentVelocity * _gameManager.TimeScale / _gameManager.SpaceScaleKm; }
        set { _currentVelocity = value * _gameManager.SpaceScaleKm / _gameManager.TimeScale; }
    }

    public bool IsVelocityIncreasing { get; private set; } = true;
    #endregion

    void Awake()
    {
        var spaceObjects = GameObject.FindGameObjectsWithTag("SpaceObject").Where(o => o.gameObject != gameObject).ToArray();
        _rigidbody = GetComponent<Rigidbody>();

        InitializeGravity(spaceObjects, _mass);
    }

    void Update()
    {
        //This code would go on a game controller and affect desired objects.
        if (Input.GetKeyDown(KeyCode.Space))
        {
            if (!_disableGravityAffected)
            {
                SetOrbit();
            }
        }
    }

    void FixedUpdate()
    {
        try
        {
            if (!initializationComplete)
            {
                Initialization();
            }

            if (!_disableGravityAffected)
            {
                _rigidbody.velocity += GetExternalNetGravityVector();
            }

            UpdateVelocityStatus();
        }
        catch (Exception e)
        {
            throw new SystemException("Problem obtaining gravity: " + e);
        }
    }

    private void UpdateVelocityStatus()
    {
        if (CurrentVelocity != _rigidbody.velocity.magnitude)
        {
            IsVelocityIncreasing = (CurrentVelocity < _rigidbody.velocity.magnitude) ? true : false;
        }

        CurrentVelocity = _rigidbody.velocity.magnitude;
    }

    private void Initialization()
    {
        _initialVelocity = GetInitialVelocity(_velocityType);
        _rigidbody.velocity = _initialDirection.normalized * (float)_initialVelocity;

        _initialVelocity = _initialVelocity * _gameManager.SpaceScaleKm / _gameManager.TimeScale;

        initializationComplete = true;
    }


    private double GetInitialVelocity(InitialVelocity initialVelocityType)
    {
        if (_spaceObjects.Length == 0)
        {
            throw new SystemException("No SpaceObjects found. Need at least 1 to get velocity");
        }

        _orbitAround = _spaceObjects
          .OrderByDescending(o => o.GetComponent<SpaceObject>().GetGravitationalPullForce(
          (o.transform.position - transform.position).magnitude))
          .First();

        var distance = (_orbitAround.transform.position - transform.position).magnitude;
        var velocity = 0.0;

        switch (initialVelocityType)
        {
            case InitialVelocity.CircularOrbit:
                velocity = _orbitAround.GetComponent<SpaceObject>().GetVelocityForCircularOrbit(distance) + _orbitAround.GetComponent<Rigidbody>().velocity.magnitude;
                break;
            case InitialVelocity.EscapeOrbit:
                velocity = _orbitAround.GetComponent<SpaceObject>().GetEscapeVelocity(distance);
                break;
            default: //Free velocity                
                velocity = _initialVelocity * _gameManager.TimeScale / _gameManager.SpaceScaleKm;
                break;
        }

        return velocity;
    }

    private void SetOrbit()
    {
        var target = _spaceObjects.First().GetComponent<SpaceObject>();
        target.GetOrbitFromCurrentVelocity(gameObject, ref _orbit);
    }
}

The code added here only affects the user input to create the orbit. This could easily be in some type of game controller that only affects the desired objects. How this is triggered would be completely project-depended. I decided to keep it here for the sake of simplicity, though it should be somewhere else. As this is global to every object with an Orbit, triggering it will display the orbits of every single object.

Unity setup

The setup to have the orbit working is very simple. We need a child Orbit GameObject to the orbiting object
Add Orbit to SpaceObject script

The orbit object just needs 2 components: A line renderer and the Orbit script:

Components required for the Orbit GameObject


Now that the setup is finished let's view a few GIFs with elliptical and hyperbolic trajectory predictions. The blue line is the object's trail. The brown line is the predicted trajectory.

Circular trajectory prediction

Elliptic trajectory prediction

Hyperbolic trajectory prediction

All these examples are in the XZ plane, to keep it simple, but they work in any other plane.

Conclusions

This is an upgrade from the previous blog entry to add One-Body problem solutions for any single object. The implementation for N-Body problems is more complicated and maybe one day I will look into it. 

As with the previous post a sample project can be found on this download link with all the relevant code.
Unity version: 2020.2.7f1

The next post will look into painting a force map using CPU, Computational and Vertex shader.

References  

[1] Keeton, C. (2014). Principles of Astrophysics Using Gravity and Stellar Physics to Explore the Cosmos. https://doi.org/10.1007/978-1-4614-9236-8

[2] Walter, U. (2018). Astronautics The Physics of Space Flight. https://doi.org/10.1007/978-3-319-74373-8

Images:

[3] Robert A. Braeunig, Figure 4.1, digital image, accessed 16 June 2021, < http://www.braeunig.us/space/orbmech.htm>

[4] Angular Parameters of Elliptical Orbit, 2005, digital image, accessed 15 June 2021, < https://en.wikipedia.org/wiki/N-body_problem#/media/File:Angular_Parameters_of_Elliptical_Orbit.png>
By at No comments:
Labels: , , , ,
Location: Badalona, Barcelona, Spain

Saturday, June 12, 2021

Gravity Part 1: An update

Introduction

Artistic representation of Earth's gravitation field (Image credit: Shutterstock)

It has been a while since I last added a new post entry. For the past year (pandemic at full) I did lots of testing regarding my physics-based gravity project.

This entry is the first of a series of 3 in which I will update what have been done so far:

The update

The code described here has basically the same functionality as the original Physics-based Solar System Lessons post. The foundations, insights and issues are described in detail in the mentioned blog entry. So, I highly recommend checking it out for references.

I did not want to change that post so I decided to create a new one with cleaner code, explaining the major changes and calculations.

I did a major code refactoring to have everything clearer and more cohesive. The code is divided into 2 blocks:
  • Generic:
    • GameManger.cs: Controls the time and space scale logic.
    • Constants.cs: Controls generic constants/ enums used in the code.
  • Gravity implementation:
    • SpaceObject.cs: This is the core logic used by the objects.
    • Gravity.cs This is the core to calculate gravity.

Generic:

The GameManager.cs class has changed quite a bit on content to only use what is important:
  • To modify the fixed delta time, everything else was not necessary and generated issues when used.
  • To calculate a “new” gravitational constant. Previously I did all the force calculations using the original units N*m^2/kg^2. Now, the distance and time is calculated using the scaled units defined in GameManager.cs. The advantage is that we would never run out of precision, regardless of the actual values (just imagine if we want to use this to simulate gravity movement in the Milky Way).
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
using System;
using UnityEngine;

public class GameManager : MonoBehaviour
{
    [Header("Scaling options")]
    [SerializeField]
    [Tooltip("1 unit = SpaceScale kilometers")]
    private float _spaceScale = 6000;

    [SerializeField]
    [Tooltip("1 unit = TimeScale seconds")]
    private float _timeScale = 100000; //1 unity unit = timeScale seconds

    [SerializeField]
    [Tooltip("Fixed update loop time (default = 0.02)")]
    private float _modifiedFixedDeltaTime = 0.001f;

    public float SpaceScaleKm => _spaceScale;
    public float TimeScale => _timeScale;

    public double GravitationalConstantUnityScaled { get; private set; }
    public float ModifiedDeltaTime => _modifiedFixedDeltaTime;

    private void Awake()
    {
        Time.fixedDeltaTime = ModifiedDeltaTime;
        GravitationalConstantUnityScaled = Constants.GRAVITATIONAL_CONSTANT_KM * Math.Pow(TimeScale, 2) / Math.Pow(SpaceScaleKm, 3);
    }
}

The Constant.cs class has increased to add generic enums used throughout the code. 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
public static class Constants
{
    public static readonly double GRAVITATIONAL_CONSTANT_KM = 6.67430e-20; //N*km^2/kg^2
    public static readonly double EPSILON = 0.01;

    public const int FREE_VELOCITY = 0;
    public const int CIRCULAR_ORBIT_VELOCITY = 1;
    public const int ESCAPE_ORBIT_VELOCITY = 2;
	
    public enum InitialVelocity
    {
        Free = FREE_VELOCITY,
        CircularOrbit = CIRCULAR_ORBIT_VELOCITY,
        EscapeOrbit = ESCAPE_ORBIT_VELOCITY
    }
}

Gravity Implementation:

Originally, I had 3 different classes (GravitationalForces, SpaceObject and Satellite). Of these, “ Satellite” is a specific case of SpaceObject and the logic between the other two is entangled. So now everything has been refactored into 2 classes:
  • Gravity.cs: The logic that does all the gravity calculation.
  • SpaceObject.cs: Any object that affects or is affected by gravitation.
For an object to be affected by this gravity implementation I am using the tag: SpaceObject. I am using the same name as the script name so it is easy to connect both.

SpaceObject Tag
Tag that enables the object to be affected by gravity

The implementation of Gravity.cs class now only takes care of functions to obtain gravitation pull force and velocities from gravitational situations.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
using System;
using UnityEngine;

public class Gravity : MonoBehaviour
{
    protected GameObject[] _spaceObjects;    
    protected GameManager _gameManager;
    private double _massSource;

    /// <summary>
    /// Get gravitational pull force from this object, depending on the distance from it
    /// </summary>
    /// <param name="distanceFromObject"> distance between 2 objects in unity units </param>
    /// <returns></returns>
    public double GetGravitationalPullForce(double distanceFromObject)
    {
        return _gameManager.GravitationalConstantUnityScaled * _massSource / Math.Pow(distanceFromObject, 2);
    }

    /// <summary>
    /// Get the needed velocity to have a stable circular orbit around a gravitational object
    /// </summary>
    /// <param name="distanceFromObject"> distance between 2 objects in unity units </param>
    /// <returns></returns>
    public double GetVelocityForCircularOrbit(double distanceFromObject)
    {
        var velocity = Math.Sqrt(_gameManager.GravitationalConstantUnityScaled * _massSource / distanceFromObject);
        return velocity;
    }

    /// <summary>
    /// Get the needed velocity to escape a gravitational object's gravity
    /// </summary>
    /// <param name="distanceFromObject"> distance between 2 objects in unity units </param>
    /// <returns></returns>
    public double GetEscapeVelocity(double distanceFromObject)
    {
        return Math.Sqrt(2 * _gameManager.GravitationalConstantUnityScaled * _massSource / distanceFromObject);
    }

    protected void InitializeGravity(GameObject[] spaceObjects, double mass)
    {
        _spaceObjects = spaceObjects;
        _massSource = mass;

        _gameManager = GameObject.Find("GameManager").GetComponent<GameManager>();
    }

    protected Vector3 GetExternalNetGravityVector()
    {
        var netForce = new Vector3();

        foreach (var spaceObjects in _spaceObjects)
        {
            netForce += GetGravity(spaceObjects) * Time.fixedDeltaTime;
        }

        return netForce;
    }

    private Vector3 GetGravity(GameObject spaceObjects)
    {
        var direction = spaceObjects.transform.position - gameObject.transform.position;
        var gravity = spaceObjects.GetComponent<SpaceObject>().GetGravitationalPullForce(direction.magnitude);

        return direction.normalized * (float)gravity;
    }
}
  
The class SpaceObject.cs on the other hand takes care of the initial changes to the object (initial velocities) and updating the forces applied to the object. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
using System;
using System.Linq;
using UnityEngine;
using static Constants;

public class SpaceObject : Gravity
{
    #region variables
    [Header("Velocity settings")]
    [SerializeField] private Vector3 _initialDirection = new Vector3(0, 0, 1);
    [SerializeField] private InitialVelocity _velocityType = InitialVelocity.Free;
    [SerializeField] private double _initialVelocity = 0; // In standard km/s
    
    [Header("SpaceObject Properties")]

    [SerializeField]
    [Tooltip("external gravity forces do not affect this object")]
    private bool _disableGravityAffected = false;

    [SerializeField]
    [Tooltip("Kilograms")]
    private double _mass = 0; //Kg

    [SerializeField]
    [Tooltip("Kilometers")]
    private double _radius = 0; //km

    [SerializeField]
    [Tooltip("Kilometers. Distance at which entry to atmosphere happens")]
    private double _entryInterfaceDistance = 0; //Km. distance at which entry to atmosphere happens

    private Rigidbody _rigidbody;
    private GameObject _orbitAround;
    private bool initializationComplete = false;

    [Header("Debug")]
    [SerializeField] private double _currentVelocity = 0; // In standard km/s
    #endregion

    #region Properties
    public double Mass => _mass;

    public double CollisionDistance { get => (_radius + _entryInterfaceDistance) / _gameManager.SpaceScaleKm; }

    public double CurrentVelocity
    {
        get { return _currentVelocity * _gameManager.TimeScale / _gameManager.SpaceScaleKm; }
        set { _currentVelocity = value * _gameManager.SpaceScaleKm / _gameManager.TimeScale; }
    }

    public bool IsVelocityIncreasing { get; private set; } = true;
    #endregion

    void Awake()
    {
        var spaceObjects = GameObject.FindGameObjectsWithTag("SpaceObject").Where(o => o.gameObject != gameObject).ToArray();
        _rigidbody = GetComponent<Rigidbody>();

        InitializeGravity(spaceObjects, _mass);
    }

    void Update()
    {

    }

    void FixedUpdate()
    {
        try
        {
            if (!initializationComplete)
            {
                Initialization();
            }

            if (!_disableGravityAffected)
            {
                _rigidbody.velocity += GetExternalNetGravityVector();
            }

            UpdateVelocityStatus();
        }
        catch (Exception e)
        {
            throw new SystemException("Problem obtaining gravity: " + e);
        }
    }


    private void UpdateVelocityStatus()
    {
        if (CurrentVelocity != _rigidbody.velocity.magnitude)
        {
            IsVelocityIncreasing = (CurrentVelocity < _rigidbody.velocity.magnitude) ? true : false;
        }

        CurrentVelocity = _rigidbody.velocity.magnitude;
    }

    private void Initialization()
    {
        _initialVelocity = GetInitialVelocity(_velocityType);
        _rigidbody.velocity = _initialDirection.normalized * (float)_initialVelocity;

        _initialVelocity = _initialVelocity * _gameManager.SpaceScaleKm / _gameManager.TimeScale;

        initializationComplete = true;
    }


    private double GetInitialVelocity(InitialVelocity initialVelocityType)
    {
        if (_spaceObjects.Length == 0)
        {
            throw new SystemException("No SpaceObjects found. Need at least 1 to get velocity");
        }

        _orbitAround = _spaceObjects
          .OrderByDescending(o => o.GetComponent<SpaceObject>().GetGravitationalPullForce(
          (o.transform.position - transform.position).magnitude))
          .First();

        var distance = (_orbitAround.transform.position - transform.position).magnitude;
        var velocity = 0.0;

        switch (initialVelocityType)
        {
            case InitialVelocity.CircularOrbit:
                velocity = _orbitAround.GetComponent<SpaceObject>().GetVelocityForCircularOrbit(distance) + _orbitAround.GetComponent<Rigidbody>().velocity.magnitude;
                break;
            case InitialVelocity.EscapeOrbit:
                velocity = _orbitAround.GetComponent<SpaceObject>().GetEscapeVelocity(distance);
                break;
            default: //Free velocity                
                velocity = _initialVelocity * _gameManager.TimeScale / _gameManager.SpaceScaleKm;
                break;
        }

        return velocity;
    }
}

These 2 classes have taken a complete refactoring to have a better segregation at what each class is supposed to do.

Sample project

A sample project can be found on this download link with all the relevant code. 
Unity version: 2020.2.7f1
By at No comments:
Labels: , , , ,
Location: Badalona, Barcelona, Spain
Subscribe to: Posts (Atom)

Gravity Part 3: Create a force map and display it (CPU, Compute Shader and Fragment Shader)

Introduction This is the third and last post in the series I am creating, based on the logic from Physics-based Solar System Lessons : ...