¿Se supone que debo modificar la constante gravitatoria con escala y por qué los cambios de fps y escala de tiempo hacen que mi órbita se rompa?

Estoy usando Unity3D y escribiendo C#

• Masa de la Tierra: 5.972e+22
• Velocidad de la Tierra: [0,0,0]
• Radio de la Tierra: 63.71
• Masa ISS: 4194.55
• Velocidad inicial de la ISS: [76.50, 0, 0]
• Altitud ISS: 4
• Constante gravitacional: 6.67408e-15

Todos los valores son los valores de la vida real / 100, ya que estoy reduciendo las cosas. (por ejemplo, la constante gravitacional fue 6.67408e-11 / 11 = 6.67408e-15.

Aquí está el código...

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;
}

Para que esto funcione, necesito establecer la constante gravitacional en: 1.59508e-21

Aquí hay una imagen de la órbita, la línea rosa es pasada, no futura (los elementos orbitales calculados que se muestran probablemente sean incorrectos):

Aquí hay una imagen de la órbita después de acelerar el tiempo:

Edición 1:

Usando esta fórmula para la aceleración y la primera parte de la respuesta de Russell Borogove, pude lograr una órbita circular usando estos valores:

• Masa de la Tierra: 5.972e+22
• Velocidad de la Tierra: [0,0,0]
• Radio de la Tierra: 63.71
• Masa ISS: 4194.55
• Velocidad inicial de la ISS: [0.0765, 0, 0]
• Altitud ISS: 4
• Constante gravitacional: 6.67408e-24 (desplazada -15 y luego +2 órdenes de magnitud)

y este código:

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;
}

Edición 2: Se agregaron cálculos de gravedad usando el script del acumulador de Russell Borogove a continuación...

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;
    }
}

Edición 3: Script de elementos orbitales actualizado (probablemente no 100% correcto)

    //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;

Edición más reciente: se solucionó el problema que tenía con la excentricidad incorrecta: estaba calculando incorrectamente la posición relativa.