Space Based Physics

Question

How would I go about about adding space-based physics like orbital mechanics into a game engine like Unreal Engine 4 or Unity 5? The goal is to generate celestial bodies procedurally: first black holes, start around those black holes, then planetary systems around those stars, then in every instance calculate the needed velocity to make stable orbits(the system would calculate the best orbit and have cases for star systems with multiple suns) them data in the format time: position will be stored in an array for each celestial body and the body will then be rendered there when needed.

Answer 1

I'm actually working on a 2D project similar to what your doing, definitely less complex in scope though.

I am also a beginner at this but can give you what I've learned so far.

You're going to need to figure out your gravity mechanics first, universe generation is easier and will have to come second.

Since I'm assuming you're making a game and not a simulation, you're going to probably use two-body gravity not n-body, the more bodies effecting each other at once the more processing time youre going to be eating. You're going to have to decide how complex you want it.

If you use Unity, you're going to want to put your physics updates in

FixedUpdate()

and you will have to write your own gravity. Using Unity's rigidbody and rigidbody.addforce() will lead to problems since unity caps rigidbody velocity.

You're going to need to scale things down a lot too.

It will also probably be best to keep your black holes, Suns, planets, and moons on the rails, once you have an orbital path calculated once just move it along that path without calculating gravity anymore. (This lets you draw the orbital path without running math constantly)

I'm on mobile at the moment but when I get a chance I'll edit this to include some links to resources and some code I've done to go from the 2 state vectors to the 6 orbital elements, but hopefully someone else will be able to help you out more, I suggest you keep the space exploration stack exchange in mind as well.

Edit 1... Here are some very useful PDFs that show you how to calculate the orbital elements, then convert them back at a specific t time. I have only used the first one in my project.

State vectors to orbital elements

Orbital elements back to state vectors at given time

Edit 2... (I won't have my code until I get home) For your forces, you're going to need torque, and velocity, I'm still using Unity's addTorque() because I'm in 2D but you may want your own. As for velocity, do not use Unity's addForce(). Store velocity as a vector3 and move the position of your object by that vector each time step. You can add any velocity changes to that vector. I'll be able to get some code on gravity and the orbital elements this evening.

You're going to need to pick a scale and stick to it, as seen in my questions on space exploration losing track of the units and not converting values properly leads to problems.

Edit 3...

Here is my code for gravity:

public Vector3 Gravity(Vector3 pos1, float m1, Vector3 pos2, float m2) {
    //Relative distance vector 
    Vector3 dvec = (pos2 - pos1);

    //Relative distance
    float d = dvec.magnitude;

    //Direction of force
    Vector3 dnorm = dvec.normalized;

    //Distance Squared
    float dsqrd = Mathf.Pow(d, 2);

    //Mass of both bodies
    float M = m1 + m2;

    //Gavitational pull
    float g = World.G * M / dsqrd;

    //Acceleration due to gravity vector
    Vector3 f = dnorm * g;

    return f;
}

Edit 4...

And finally, here is my huge script that finds a bunch of sciency stuff for gravity. Some of it you will ignore, I can't confirm all of it is 100% correct, but you can use it as a resource as I suffered through it myself and wouldn't wish the same torture on my worst enemy

public OrbitData CalculateOrbitalData(Vector3 pos1, Vector3 vel1, float m1, Body body) {
    if (body == null) {
        return new OrbitData();
    }

    Vector3 pos2 = body.position;
    Vector3 vel2 = body.velocity;
    float m2 = body.mass;

    //Relative Position Vector
    Vector3 r = pos1 - pos2;

    //Distance between bodies
    float rmag = r.magnitude;

    //Relative Velocity Vector
    Vector3 v = vel1 - vel2;

    //Velocity Magnitude
    float v1 = vel1.magnitude;

    //Relative Velocity Magnitude
    float vm = v.magnitude;

    //Specific Angular Momentum
    Vector3 h = Vector3.Cross(r, v);

    //Standard Gravitational Parameter
    float µ = World.G * m2;

    //Eccentricity Vector
    Vector3 evec = (Vector3.Cross(v, h) / µ) - (r / Vector3.Magnitude(r));

    //Eccentricity
    float e = Vector3.Magnitude(evec);

    //Vector to Ascending Node
    Vector3 n = new Vector3(-h.x, h.z, 0);

    //True Anomaly
    float ta = Mathf.Acos((Vector3.Dot(evec, r)) / (e * Vector3.Magnitude(r)));
    if (Vector3.Dot(r, v) < 0)
        ta = (2 * Mathf.PI) - ta;

    //Longitude of Ascending Node (2D)
    float Ω = 0;

    //Inclination (2D)
    float i = 0;

    //Argument of Periapsis
    float ω = Mathf.Atan2(evec.y, evec.x);
    float ωdegrees = ω * (180 / Mathf.PI);
    if (e == 0) {
        ω = 0;
        ωdegrees = 0;
    }

    //Eccentric Anomaly
    //float E = 2 * Mathf.Atan2(Mathf.Tan(ta / 2), Mathf.Sqrt((1 + e) / (1 - e)));
    float E = 2 * Mathf.Atan(Mathf.Tan(ta / 2) / Mathf.Sqrt((1 + e) / (1 - e)));

    //Mean Anomaly 
    float M = E - (e * Mathf.Sin(E));

    //Semi-Major Axis
    float a = 1 / ((2 / Vector3.Magnitude(r)) - (Mathf.Pow(Vector3.Magnitude(v), 2) / µ));

    //Apoapsis
    float ap = a * (1 + e);

    //Periapsis
    float pe = a * (1 - e);

    //Orbital Period
    float T = (2 * Mathf.PI) * Mathf.Sqrt(Mathf.Pow(a, 3) / µ) / 60;
    float fT = T * 1.205f;

    //Mean motion
    float m = (2 * Mathf.PI) / fT;

    //Time since apoapsis
    float tsa = Mathf.Sqrt(Mathf.Pow(a, 3) / µ) * ((E + Mathf.PI) - e * Mathf.Sin((E + Mathf.PI))) / 60 * 1.205f;

    //Time to apoapsis
    float tta = fT - tsa;

    //Time since periapsis
    float tsp = tsa - fT / 2;
    if (tsp < 0)
        tsp = fT / 2 + tsa;

    //Time to periapsis
    float ttp = tta - fT / 2;
    if (ttp < 0)
        ttp = fT / 2 + tta;

    //Perifocal distance
    float rp = (Vector3.Dot(h, h) / µ) / e + 1;

    //Prograde rotation angle
    float pr = 0;
    if (vel1 != Vector3.zero) pr = Quaternion.AngleAxis(Mathf.Atan2(vel1.y, vel1.x) * Mathf.Rad2Deg, Vector3.forward).eulerAngles.z - 90;

    //Retrograde rotation angle
    float rr = pr + 180;

    //Counterclockwise forward rotation angle
    float ccwR = 0;
    float rad = 0;
    Vector2 offset = transform.position - body.transform.position;
    if (offset.x != 0) rad = (Mathf.Atan(offset.y / offset.x));
    if (offset.x < 0) rad += Mathf.PI;
    ccwR = rad * (180 / Mathf.PI);

    //Clockwise forward rotation angle
    float cwR = ccwR + 180;

    //Up rotation angle
    float uR = cwR + 90;

    //Down rotation angle
    float dR = cwR - 90;

    //circularize dV 
    float circ_dV = Mathf.Sqrt(µ / ap) - Mathf.Sqrt((pe * µ) / (ap * (pe + ap) / 2));

    //circularize burn tick count
    float circ_btt = circ_dV / forwardThrust;

    //Circularize burn time
    float circ_bt = circ_btt / (1 / World.timeStep);

    return new OrbitData(m1, pos1, v1, vel1, m2, vel2.magnitude, vel2, pos2, World.G,
        r, v, rmag, h, µ, evec, e, n, ta, Ω, i, ω, ωdegrees, E, M, a, ap, tsa, tta,
        pe, tsp, ttp, orbit.forcedApoapsisPosition, orbit.forcedPeriapsisPosition,
        T, fT, rp, pr, rr, cwR, ccwR, uR, dR, circ_dV, circ_btt, circ_bt);
}