Am I supposed to modify the gravitational constant with scale and why do fps & time scale changes cause my orbit to break?

Question

I am using Unity3D and writing C#

• Earth Mass: 5.972e+22
• Earth Velocity: [0,0,0]
• Earth Radius: 63.71
• ISS Mass: 4194.55
• Starting ISS Velocity: [76.50, 0, 0]
• ISS Altitude: 4
• Gravitational Constant: 6.67408e-15

All values are the real life values / 100, as I am scaling things down. (ex. Gravitational Constant was 6.67408e-11 / 11 = 6.67408e-15.

Here is the code...

void Gravity(Gravity body) {
    float dSquared = (body.transform.position - transform.position).sqrMagnitude;
    float force = (body.G * shipMass * body.planetMass) / dSquared;
    Vector3 forceDirection = (body.transform.position - transform.position).normalized;
    Vector3 forceVector = (forceDirection * force);
    shipVelocity += forceVector * Time.deltaTime;

    transform.position += shipVelocity * Time.deltaTime;
}

In order to make this work I need to set the gravitational constant to: 1.59508e-21

Here is an image of the orbit, the pink line is past not future path (the calculated orbital elements shown are most likely wrong):

Here is an image of the orbit after I speed up time:

Edit 1:

Using this formula for acceleration and the first part of Russell Borogove's answer I was able to achieve a circular orbit using these values:

• Earth Mass: 5.972e+22
• Earth Velocity: [0,0,0]
• Earth Radius: 63.71
• ISS Mass: 4194.55
• Starting ISS Velocity: [0.0765, 0, 0]
• ISS Altitude: 4
• Gravitational Constant: 6.67408e-24 (Shifted -15 then +2 orders of magnitude)

and this code:

void Gravity(Gravity body) {
    float dSquared = (body.transform.position - transform.position).sqrMagnitude;
    float M = body.planetMass;
    Vector3 forceDirection = (body.transform.position - transform.position).normalized;
    float g = G * M / dSquared;
    Vector3 forceVector = (forceDirection * g);
    shipVelocity += forceVector * Time.deltaTime;
}

Edit 2: Added gravity calculations using Russell Borogove's accumulator script below...

void Update() {
    Gravity(orbitee);
    ApplyForce();
}

void ApplyForce() {
    transform.position += shipVelocity * Time.deltaTime;
    transform.Rotate(transform.forward * shipTorque * Time.deltaTime);
}

void Gravity(Gravity body) {
    accumulator += Time.deltaTime;
    while (accumulator > 0.0f) {
        float dSquared = (body.transform.position - transform.position).sqrMagnitude;
        float M = body.planetMass;
        Vector3 forceDirection = (body.transform.position - transform.position).normalized;
        float g = G * M / dSquared;
        Vector3 forceVector = (forceDirection * g);
        shipVelocity += forceVector * simtime;
        accumulator -= simtime;
    }
}

Edit 3: Updated Orbital Elements Script (Probably not 100% correct)

    //Mass of satellite
    float m1 = shipMass;

    //Mass of planet
    float m2 = orbiting.planetMass;
    float M = m2;

    //Relative Position Vector
    Vector3 r = FindRelativePosition(orbitee.transform.position, transform);

    //Relative Velocity Vector
    Vector3 v = shipCurVel;
    //#!# add relative to planet velocity

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

    //Standard Gravitational Parameter
    float µ = G * (m1 + m2);
    //float µ = G * M;

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

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

    //Vector to Ascending Node #!#
    Vector3 n = Vector3.Cross(new Vector3(0, 0, 1), h);

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

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

    //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);
    //float ω = Mathf.Acos((Vector3.Dot(n, evec)) / (e * Vector3.Magnitude(n)));

    //Mean Anomaly 
    float MA = 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;

Most recent edit: fixed issue I was having with eccentricity being incorrect - was calculating relative position wrong.

Answer 1

Big G is 6.67408e-11 m³ kg^-1 sec^-2.

So if you uniformly use 100km units instead of m units, that's 5 magnitudes different, so G has to change by 15 magnitudes because the linear dimension is cubed.

Mass units are 100kg, G uses -1 of those, so G goes the other way by 2 magnitudes.

I have no idea what you're doing with seconds. If seconds are seconds (and it takes 90 minutes to complete an orbit), then there's no change to G. If you're using some other conversion factor, then double the number of zeros you're shifting seconds by to change G.

As for the instability, that's inevitable when you do time-step simulations of continuous functions. You can split your Unity Update time slice into some small fixed quantum and run your simulation function several times for each game frame:

float accumulator = 0.0f;
float simtime = 0.001f;

void Update( void )
{
    accumulator += Time.deltaTime;
    while (accumulator > 0.0f)
    {
         SimulateGravityForTime( simtime );
         accumulator -= simtime;
    }
}

This runs your sim once for each millisecond of time elapsed in the frame. That should keep things pretty stable. The accumulator is there to retain fractional milliseconds from frame to frame. You'll spend more CPU time doing repeated computations, but a modern computer should be able to afford it easily.

---

Regarding the units, if I were in your shoes I would separate the view (based on Unity transforms) from the model (the physics simulation). This would allow you to run the simulation entirely in standard units, eliminating much confusion and reducing the likelihood of errors. Inside the physics model you'd always use m/kg/s units and keep track of positions in vectors that are not part of Unity transforms. Once per frame you'd carefully convert to Unity units, and set the positions of the transforms.