2
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ You're using ship mass to compute force, but not dividing by mass when you apply it as acceleration (a quick check, 1.595 (mantissa of your G) x 4194 (mass of ship) yields about 6700 (matching mantissa of correct G) so that's likely the problem). It's generally easier to skip computing force and just compute acceleration, which is independent of ship mass. And again you're mixing up 100m and 100km units. $\endgroup$ Jun 16, 2016 at 3:30
  • $\begingroup$ I will give that a try and report back. So the velocity should be 0.0765 instead of 76.5? $\endgroup$
    – Wafer
    Jun 16, 2016 at 3:37
  • $\begingroup$ It's probably okay if you get the right number of zeros in your G. $\endgroup$ Jun 16, 2016 at 3:40
  • $\begingroup$ Does G change with scale or is it something that should be constant? I had to make it significantly smaller compared to the change in all the other parameters. $\endgroup$
    – Wafer
    Jun 16, 2016 at 3:46
  • $\begingroup$ Of course it does. Elaborated in answer. $\endgroup$ Jun 16, 2016 at 3:57

1 Answer 1

7
$\begingroup$

Big G is 6.67408e-11 m3 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.

$\endgroup$
11
  • 1
    $\begingroup$ +1 for great answer, +1 for including a little code, -1 for aiding and abetting creative and individualized physical units without even a little caveat about what can go wrong down the road. :) Even inches vs centimeters almost killed Hubble. $\endgroup$
    – uhoh
    Jun 16, 2016 at 5:27
  • 1
    $\begingroup$ I think I've been as clear in the discussions I've had with OP (in comments in this and the previous question) that paying attention to the units is important. That said, I'll add a bit about view-model separation. $\endgroup$ Jun 16, 2016 at 13:20
  • 1
    $\begingroup$ Good - I'll start working on a way to change my vote to +2 which your answer deserves. It involves integer overflow and something about row hammering and then some tachyons, but I'm just not there yet. (humor) $\endgroup$
    – uhoh
    Jun 16, 2016 at 13:28
  • $\begingroup$ @RussellBorogove I've finished the gravity and used your accumulator script here. When I convert T to minutes you can see that it matches the 90m orbital period here. I have edited my post to include the orbital elements script. Does there appear to be anything wrong with my orbital elements? I'm sure the AP and PE are incomplete or wrong as I need to modify them by radius to get numbers that are close, or the magnitude of r in order to get accurate numbers at t: 0. $\endgroup$
    – Wafer
    Jun 17, 2016 at 3:01
  • 1
    $\begingroup$ The 2 in Atan2 in computer math libraries refers to two arguments, not multiplying by 2. One-argument Atan can't tell the difference between ( + / + ) and ( - / - ) inputs, which is important in some trig applications -- I don't know if this is one of them. The second of your options looks right, the first needs a 2 *. If your results are sane for half the orbit and insane for the other half, you need Atan2. $\endgroup$ Jun 17, 2016 at 4:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.