Physics-based Solar System Lessons

I have always been a fan of space simulators and physics. I thought about the possibility of creating a physics-based space simulator. I decided to create a simple prototype using Newton’s gravitational formula to move objects in space. The idea was to start with the major planets and see how they interacted with each other, as well as with other objects (asteroids, spaceships, missiles, etc.).

Results

As the original idea was a bit too complex I based my prototype on visualizing Earth-Moon system (and a couple of satellites) and the inner solar system (no satellites). From here everything can be as complex as desired.

Here are the results of Earth-Moon system with 2 satellites orbiting (one orbits Earth the other the Moon):

Distance scale = 6000 km, Time scale = 500000 seconds

I added some trail Renderer to visualize the movement of the objects in space.

Additionally I also created the inner solar system (Mercury, Venus, Earth and Mars)

Distance scale = 1e+07 km, Time scale = 3200000 seconds

Even though these are top view demonstrations the implementation is done in 3D, allowing for more interesting camera views or orbits inclinations. These are very nice results, but there were many things learnt and discovered.

The Math

In order to create this physically realistic system we need to use Newton's formulas. Its gravitation formula works great and does not contain any complicated math.

To see how much force is applied to an object at any given position (regardless of mass) you just apply Newton's law of universal gravitation and Newton’s second law together.

Newton's law of universal gravitation:

\[F = G\dfrac{m_1 m_2}{r^2}\]

Where:

G = Gravitational Constant = 6.67430e-11 N*m^2/kg^2
m_1 = Mass Object_1 in kg (e.g. Satellite)
m_2 = Mass Object_2 in kg (e.g. Earth)
r = distance from Object_1 to Object_2 in meters

Newton’s second law:

\[F = m_1 a\]

Where

a = acceleration in m/s^2

This will provide the acceleration, which is what we need to apply to our object (e.g. Satellite):

\[a = G \dfrac{m_2}{r^2}\]

The problems

I soon realized that there were a few things that were giving me issues:

Precision issues (Size of the Solar System).
Prediction of orbits.
Speed up planetary movement.

After doing some research I found a video from Unite 2013 Kerbal Space Program that ran into most of these issues. Recommended view for anyone interested in this kind of simulation.

Precision issues

Unity works with Vector2/ Vector3 and these components are floats. Size of solar system is huge. Float only allows 7 digits of accuracy. Details on this can be read in Unity - Coordinates and scales.

Now use Pluto’s orbit as reference example. Pluto’s aphelion distance is 7.37593e+12. If you plug this value in the position of an object you will see from the message given by Unity that there might be issues with these kinds of large values:

This message appears when units are over 100000, or 1e+5

Solution: apply distance scaling, e.g. 1 unity unit = 6000 km.

Prediction of orbits

Considering multiple massive objects, applying newton’s gravitational formula made prediction complicated as forces affected the orbit differently depending on position and time. This is relevant when an object gets close enough to a massive source. A good example of orbit modification would be this Earth-Moon-Satellite example:

Moon the blue orbit, satellite the brown orbit

Solution: Only consider the most massive object affecting (usually the sun in solar scale or a planet if close to one).

Speed up planetary movement

To visualize orbits in real time is simple, but extremely slow. I.e. The moon takes almost a month to orbit the earth, would anyone wait for a month to view the actual orbit?

One can scale up unity timespan up to 100, which is not enough. Optionally one can scale up accelerations and velocities with a timestep to simulate increased speed. This generates noise in the actual orbit. See the 2 gifs for comparison.

Time scale = 5000. Precise.

If we speed up the time (time scale = 50000) we run into the following issue:

Time scale = 50000. Orbit not precise.

One might wonder why, but it has a very simple answer: Discretization of a continuous function distorts the signal. This means in our specific case the timestep AND refresh rate are too high and affects the precision in which the forces are calculated.

Solution 1: Do not scale more than certain amount. Not viable for large systems (e.g. Solar System)

Solution 2: Modify Time.DeltaTime to create more calls between FixedDeltaTime. This brings the issue of performance into account. So be careful not to abuse this option.

I am using a Time.DeltaTime = 0.001 and results are good.

The implementation

There are many ways to implement this and my solution is a quick prototype among many possibilities. You are welcome to use my implementation as a baseline.

The simplest setup is to have 3 different scripts:

GravitationalForces: The logic that does all the gravity calculation.
SpaceObject: This will belong to objects that create a gravitational pull on other objects.
Satellite: This would be a “weightless” type of object compared to SpaceObject. A few tones (or hundreds) vs billions of tones.

I also have 2 more classes for control:

Constants: Where I keep any constants I might be using.
GameManager: Here I control the time (included fixed delta time – explained later why) and space scale. I also control the speeding up of unity internal clock.

Code

Here are the sample code snippets:

GravitationalForces.cs

using System;
using UnityEngine;
using System.Linq;

public class GravitationalForces : MonoBehaviour
{
    protected const int FREE_VELOCITY = 0;
    protected const int CIRCULAR_ORBIT_VELOCITY = 1;
    protected const int ESCAPE_ORBIT_VELOCITY = 2;

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

    [Header("Velocity configuration")]

    [SerializeField]
    [Tooltip("Auto calculate velocity of the object: \n" +
        " - Free: No constrain. \n" +
        " - Circular: Orbit around closest massive object. \n" +
        " - Escape: Leave orbit from closest massive object.")]
    protected InitialVelocity _velocityType;

    [SerializeField]
    [Tooltip("Velocity in km/s")]
    protected double _initialVelocity = 0;//meters/second
    [SerializeField]
    protected Vector3 _initialDirection = new Vector3();

    //[Header("Debug")]
    //[SerializeField]
    protected GameObject[] _spaceObjects;
    protected Rigidbody _rigidbody;
    protected GameManager _gameManager;

    protected bool initializationComplete = false;

    protected void Init()
    {
        _rigidbody = GetComponent<Rigidbody>();
        _gameManager = GameObject.Find("GameManager").GetComponent<GameManager>();
        _spaceObjects = GameObject.FindGameObjectsWithTag("GravitationalSpaceObject").Where(o => o.gameObject != gameObject).ToArray();

        _initialVelocity = GetVelocity((int)_velocityType);
        _rigidbody.velocity = _initialDirection * (float)_initialVelocity * Constants.KM_TO_METERS * _gameManager.TimeScale / _gameManager.SpaceScaleMeters;

        initializationComplete = true;
    }

    protected void ApplyGravity()
    {
        var gravityForces = new Vector3();

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

        _rigidbody.velocity += gravityForces;
    }

    protected Vector3 GetGravity(GameObject spaceObject)
    {
        var direction = spaceObject.transform.position - transform.position;        
        var distance = direction.magnitude * _gameManager.SpaceScaleMeters;

        var gravity = spaceObject.GetComponent<SpaceObject>().GetGravitationalPullForce(distance);
        var gravityScaled = (float)(gravity * Math.Pow(_gameManager.TimeScale, 2) / _gameManager.SpaceScaleMeters);

        return direction.normalized * gravityScaled;
    }

    private double GetVelocity(int velocityType)
    {
        if (_spaceObjects.Length == 0)
        {
            throw new SystemException("No spaceObjects found. Needed at least 1 to get velocity");
        }

        var target = _spaceObjects
          .OrderByDescending(o => o.GetComponent<SpaceObject>().GetGravitationalPullForce(
          (o.transform.position - transform.position).magnitude) * _gameManager.SpaceScaleMeters)
          .First();

        var distance = (target.transform.position - transform.position).magnitude * _gameManager.SpaceScaleMeters;
        var velocity = 0.0;

        switch (velocityType)
        {
            case CIRCULAR_ORBIT_VELOCITY:
                velocity = target.GetComponent<SpaceObject>().GetVelocityForCircularOrbit(distance) / 1000 + target.GetComponent<Rigidbody>().velocity.magnitude;
                break;
            case ESCAPE_ORBIT_VELOCITY:
                velocity = target.GetComponent<SpaceObject>().GetEscapeVelocity(distance) / 1000;
                break;
            default: //Free velocity
                Debug.Log("FREE_VELOCITY (" + _velocityType + ")");
                velocity = _initialVelocity;
                break;
        }

        return velocity;
    }
}

SpaceObject.cs

using System;
using UnityEngine;

public class SpaceObject : GravitationalForces
{
    #region variables
    [Header("Space Object properties")]
    
    [SerializeField]
    private bool _onlyUseGravity = true;

    [Header("Space Object properties")]
    [SerializeField]
    private double _mass = 0;

    #endregion

    public bool OnlyUseGravity
    {
        get { return _onlyUseGravity; }
        set { _onlyUseGravity = value; }
    }

    #region Unity functions

    void FixedUpdate()
    {
        if (!_onlyUseGravity)
            return;

        if (!initializationComplete)
            Init();

        ApplyGravity();
    }

    void OnDrawGizmosSelected() //Debugging in editor
    {
        if (_spaceObjects != null)
        {
            Gizmos.color = Color.gray;

            foreach (var spaceObject in _spaceObjects)
            {
                Gizmos.DrawLine(transform.position, spaceObject.transform.position);
            }
        }
    }
    #endregion

    /// <summary>
    /// Get gravitational pull force from this object, depending on the distance from it
    /// </summary>
    /// <param name="distanceMetersFromObject"> distance between 2 objects in meters </param>
    /// <returns></returns>
    public double GetGravitationalPullForce(double distanceMetersFromObject)
    {
        return Constants.GRAVITATIONAL_CONSTANT * _mass / Math.Pow(distanceMetersFromObject, 2);
    }

    /// <summary>
    /// Get the needed velocity to have a stable circular orbit around a gravitational object
    /// </summary>
    /// <param name="distanceMetersFromObject"> distance between 2 objects in meters </param>
    /// <returns></returns>
    public double GetVelocityForCircularOrbit(double distanceMetersFromObject)
    {
        return Math.Sqrt(Constants.GRAVITATIONAL_CONSTANT * _mass / distanceMetersFromObject);
    }

    /// <summary>
    /// Get the needed velocity to escape a gravitational object's gravity
    /// </summary>
    /// <param name="distanceMetersFromObject"> distance between 2 objects in meters </param>
    /// <returns></returns>
    public double GetEscapeVelocity(double distanceMetersFromObject)
    {
        return Math.Sqrt(2 * Constants.GRAVITATIONAL_CONSTANT * _mass / distanceMetersFromObject);
    }
}

I have added a few additional functions to obtain relevant velocities more easily:

Circular orbit: It returns the needed velocity to keep a circular orbit to the Object that affects the strongest pull.

\[v=\sqrt{G\dfrac{m_1}{r}}\]

Escape from orbit: It returns the needed velocity to break free from the Object that affects the strongest pull.

\[v=\sqrt{2G\dfrac{m_1}{r}}\]

Satellite.cs

using UnityEngine;

public class Satellite : GravitationalForces
{
    // Update is called once per frame
    void FixedUpdate()
    {
        if (!initializationComplete)
            Init();

        ApplyGravity();
    }

    void OnDrawGizmosSelected()
    {
        if (_spaceObjects != null)
        {            
            Gizmos.color = Color.white;

            foreach (var spaceObject in _spaceObjects)
            {
                Gizmos.DrawLine(transform.position, spaceObject.transform.position);
            }
        }
    }
}

Constants.cs

public static class Constants
{
    public static readonly double GRAVITATIONAL_CONSTANT = 6.67430e-11; //N*m^2/kg^2
    public static readonly int KM_TO_METERS = 1000;
}

GameManager.cs

using UnityEngine;

public class GameManager : MonoBehaviour
{
    [Header("Unity Time speed (1 = realtime)")]
    [SerializeField]
    private float _timeScaleMod = 3.0f;
    [SerializeField]
    private float _timeScaleInfo;

    [Header("Scaling options")]
    [SerializeField]
    [Tooltip("1 unit = SpaceScale kilometers")]
    private float _spaceScale = 6000;

    [SerializeField]
    [Tooltip("1 unit = TimeScale seconds")]
    private float _timeScale = 100; //1 unity unit = timeScale seconds
    [SerializeField]
    [Tooltip("Fixed update loop time (default = 0.02)")]
    private float _modifiedFixedDeltaTime = 0.02f;//Default

    public float SpaceScaleMeters
    {
        get { return _spaceScale * 1000; }
    }
    public float SpaceScaleKm
    {
        get { return _spaceScale; }
    }

    public float TimeScale
    {
        get { return _timeScale; }
    }
    public float ModifiedDeltaTime
    {
        get { return _modifiedFixedDeltaTime; }
    }

    private void Awake()
    {
        Time.fixedDeltaTime = ModifiedDeltaTime;
    }

    // Update is called once per frame
    void Update()
    {
        UpdateTimeScale();
    }

    private void UpdateTimeScale()
    {

        if (Input.GetKeyDown(KeyCode.KeypadPlus))
        {
            Time.timeScale = Mathf.Clamp(Time.timeScale + _timeScaleMod, 1.0f, 99);
            _timeScaleInfo = Time.timeScale;
        }

        if (Input.GetKeyDown(KeyCode.KeypadMinus))
        {
            Time.timeScale = Mathf.Clamp(Time.timeScale - _timeScaleMod, 1.0f, 99);
            _timeScaleInfo = Time.timeScale;
        }
    }
}

Setup in Scene

For a basic setup I can show you my most basic layout. As a generic note, any object that has a rigidbody does not have gravity enabled in that component. We will simulate that.

GameManager: This is an empty GameObject where we control the distance and time scaling.

Earth: This is the center of the system and is in charge of affecting most of the gravitational pull. The Mass variable is extremely important here.

Any object with insignificant mass. It only needs to be told what type of velocity to be applied to it.

Moon: The other massive object in this case. It is important to have an initial velocity and Mass setup.

Data

One of the biggest issues was to find the proper data to obtain the gravitational pull and the orbit information. The most important thing to know is what are the positions relative from Source and orbiting objects (e.g. Sun and Earth, or Earth and Moon/Satellites).

Orbits

The orbits of objects always follow an elliptical trajectory (a circular orbit is a specific case of an ellipse).

In the case of the Sun-Earth system. An ellipse has the point where the sun is at one of the foci points (in a circular orbit that point is the center). The closest point along the major axis is the Perihelion and the farthest point is the Aphelion. If you have the velocity at either of those points and the distance you will obtain an orbit like the picture below.

From NASA website

For details on ellipse math you can check this link. For details on orbits you can check the NASA or ESA site.

For realistic representation of the solar system it is important to use orbital parameters to describe the position of the planets.

The trick is to find all the required information needed (e.g. Aphelion, Perihelion, eccentricity, inclination, etc.). That data can be found in the following links:

Space Objects fact sheets: https://nssdc.gsfc.nasa.gov/planetary/factsheet/

Mercury: https://nssdc.gsfc.nasa.gov/planetary/factsheet/mercuryfact.html

Venus: https://nssdc.gsfc.nasa.gov/planetary/factsheet/venusfact.html

Earth: https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html

Mars: https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html

Jupiter: https://nssdc.gsfc.nasa.gov/planetary/factsheet/jupiterfact.html

Saturn: https://nssdc.gsfc.nasa.gov/planetary/factsheet/saturnfact.html

Neptune: https://nssdc.gsfc.nasa.gov/planetary/factsheet/neptunefact.html

Uranus: https://nssdc.gsfc.nasa.gov/planetary/factsheet/uranusfact.html

Moon: https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html

Conclusion

After creating this prototype I have learnt that math is a beautiful tool to create unique behaviors and especially to simulate real-life situations. I have also realized there are many things that can be added to this small project. Most importantly the limitations have become easy to spot and there are solutions that can fix most of the problems.

For me the biggest issue was the time scaling math. The logic is simple but the implementation was not straightforward.

I personally would recommend working in multi-scene scenarios depending on the scale you are using. It is not the same Earth-Moon distances to Sun-Earth distances. Moons and satellites cannot be displayed in the second case.

I hope anyone who reads this will find this information useful.

Ideas - Game Development and Vision

Monday, March 16, 2020