Following in HarvesteR's footsteps (Unity maths)

[Corbald]

(Mods: Is this the right place to ask this?)

Hey all, I need some help! I have been working for a while on a project of my own. I have a lot of stuff that I have already done, but now it comes time to do the thing I've been dreading the most... the math of orbital dynamics and conics. I have a rudimentary understanding of how cones draw ellipses and how planets sit in the foci, but I'm trying to really wrap my mind around the math. I can't read mathematical notation and have only a high-school education. I'm ADHD and haven't been in school for 24 years. As you can imagine, this has been an uphill battle for me, the entire time! None the less, I have tackled a lot of difficult problems, including dynamic LoD, procedural generation, networking, optimization and more.

In trying to make the numbers actually mean something to my idiot brain, I stumbled upon a post by none other than HarvesteR, himself (http://www.orbiter-forum.com/showthread.php?t=20580). He seems to be asking some of the same questions that I am, and seems to be at around the same stage that I am (at that moment in time). Now, we all know that he was successful in his attempts, so I think I'm on the right track.

In that post, Arrowstar very kindly supplies him with a MATLAB example which spits out all the relevant data points for an orbit. I'm really hoping you guys can help me get the MATLAB gibberish to produce some real numbers, so that I can actually learn something and get my mind to visualize these concepts. I have a lot of references that I have dug up via Googling around, and I hope that if I can get this thing running, I can actually make some headway.

Here's what I have: (I have commented relevant lines of code with the pertinent questions I have.)

using UnityEngine;
using System.Collections;

public class Manager : MonoBehaviour
{

    // Use this for initialization
    void Start()
    {
        GameObject[] Objects = GameObject.FindGameObjectsWithTag("Star");

        //Initial 'kick'
        foreach (GameObject ObjectA in Objects)
        {
            ObjectA.GetComponent<Rigidbody>().AddForce(new Vector3(0, 0, 1)); //Initial 'kick'
        }
    }


    void ApplyGravity(Rigidbody A, Rigidbody 
    {
        //This is how to get the distance vector between two objects. (I totally stole most of this function, 
        //but I understand it now, so I'll rewrite it to suit my project, later.)
        Vector3 dist = B.transform.position - A.transform.position;
        float r = dist.magnitude;
        dist /= r;

        //This is the Newton's equation
        //G = 6.67 * 10^-11 N.m².kg^-2
        double G = 6.674f * (10 ^ 11);
        float force = ((float)G * A.mass * B.mass) / (r * r);

        //Then, just apply the forces
        A.AddForce(dist * force);
        B.AddForce(-dist * force);

        //note that there are only two objects in my scene, a 'planet' with x,y,z locked, and a 'moon'
        //the 'moon' is the only object which will get a force, for this example, due to the constraints on the other body.

        //DrawOrbit(A.position, A.position + A.velocity, dist * force); //this will call the orbital calculation stuff below.

    }

Here's the code for the function I'm having issues with

    void DrawOrbit(Vector3 rVect, Vector3 vVect, float muCB)
    {
        //Here begins the MATLAB code provided to HarvesteR
        //  function [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom] = getKeplerFromState(rVect,vVect,muCB) 
        //% getKeplerFromState() returns Keplerian orbital elements when provided
        //% with the state (cartesian position vector, cartesian velocity vector) of
        //% a spacecraft or celestial body.
        //%
        //% INPUTS
        //% rVect - a 3x1 vector that contains the x,y,z components of the orbiting
        //% body's current position relative to the central body.  Units: [km]
        //% vVect - a 3x1 vector that contains the x,y,z components of the orbiting
        //% body's current velocity vector relative to the central body.  Units:
        //% [km/sec]
        //% muCB - the gravitational parameter of the central body.  Units: km^3/s^2 //(What, of the above, goes here?
        //%
        //%OUTPUTS
        //% sma - semi-major axis of the orbit.  Units: [km]
        //% ecc - eccentricity of the orbit.  Units: dimensionless
        //% inc - inclination angle of the orbit.  Units: radian
        //% longAscNode - Longitude of ascending node of the orbit.  Units: radian
        //% ArgPeri - Argument of periapse of the orbit.  Units: radian.
        //% TrueAnom - Current true anomaly of the spacecraft/body in the orbit.
        //% Units: radian

        //r=norm(rVect);
        //rUnitVect=rVect/r;
        //v=norm(vVect);

        float r = rVect.magnitude; //matlab 'norm' is Unity Vector3.magnitude?
        Vector3 rUnitVect = rVect.normalized; // rVect/r is Unity Vector3.normalized?
        float v = vVect.magnitude;

        //hVect=cross(rVect,vVect);
        //h=norm(hVect);
        //hUnitVect=hVect/h;
        //ThetaUnitVect=cross(hUnitVect,rUnitVect);

        Vector3 hVect = Vector3.Cross(rVect, vVect);
        float h = hVect.magnitude;
        Vector3 hUnitVect = hVect.normalized;
        Vector3 ThetaUnitVect = Vector3.Cross(hUnitVect, rUnitVect);


        //Energy=v^2/2 - muCB/r;
        //sma=-muCB/(2*Energy);

        float Energy = Mathf.Pow( v, 2 )/ 2 - muCB / r; //Am I doing order of operation right?
        float sma = -muCB / (2 * Energy);

        //p=h^2/muCB;
        //ecc=sqrt(-p/sma + 1);

        float p = Mathf.Pow(h, 2) / muCB;
        float ecc = Mathf.Sqrt(-p / sma + 1);

        //TrueAnom=acos((p/r - 1)/(ecc));
        //if(dot(rVect,vVect)<0)
        //    TrueAnom=-TrueAnom;
        //end

        float TrueAnom = Mathf.Acos((p / r - 1) / ecc);
        if (Vector3.Dot(rVect, vVect) < 0)
        {
            TrueAnom = -TrueAnom;
        }

        //inc=acos(hUnitVect(3));

        float inc = Mathf.Acos(hUnitVect[2]); //Unity vectors are from 0-2, not 1-3

        //longAscNode_1=AngleZero2Pi(asin(hUnitVect(1)/sin(inc))); //is AngleZero2Pi a matlab function? what is this?
        //longAscNode_2=AngleZero2Pi(pi-asin(hUnitVect(1)/sin(inc)));
        //longAscNode_3=AngleZero2Pi(acos(-hUnitVect(2)/sin(inc)));
        //longAscNode_4=AngleZero2Pi(-acos(-hUnitVect(2)/sin(inc)));

        //longAscNodeSet1=round(1000*[longAscNode_1,longAscNode_2])/1000; //is this a Mathf.Round as a vector2?
        //longAscNodeSet2=round(1000*[longAscNode_3,longAscNode_4])/1000;

        //[val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2); //am I finding the intersection of two vector2's? I don't know how to read this!
        //longAscNode=longAscNodeSet1(ia); //I don't know how to read this!

        //Theta_1=AngleZero2Pi(asin(rUnitVect(3)/sin(inc))); //Quesion as above
        //Theta_2=AngleZero2Pi(pi-asin(rUnitVect(3)/sin(inc)));
        //Theta_3=AngleZero2Pi(acos(ThetaUnitVect(3)/sin(inc)));
        //Theta_4=AngleZero2Pi(-acos(ThetaUnitVect(3)/sin(inc)));

        //ThetaSet1=round(1000*[Theta_1,Theta_2])/1000; //Quesion as above
        //ThetaSet2=round(1000*[Theta_3,Theta_4])/1000;

        //[val,ia,ib]=intersect(ThetaSet1,ThetaSet2); //Quesion as above
        //Theta=ThetaSet1(ia); //Quesion as above

        //ArgPeri=Theta-TrueAnom; //If I could figure out the above, this might make sense!


        //some stuff will go here to draw the elipse/orbits. I'll tackle that later, after I make sense of the rest.
    }

Finally, some code to apply a basic orbit, it's bad, but it's working

    void FixedUpdate()
    {
        //Get every object 
        GameObject[] Objects = GameObject.FindGameObjectsWithTag("Star"); //This runs too often. Extract to Start() later to optimize.

        //the gravity between each couple of object is calculated
        foreach (GameObject ObjectA in Objects)
        {
            foreach (GameObject ObjectB in Objects)
            {
                //Objects must not self interact 
                if (ObjectA == ObjectB)
                    continue;

                ApplyGravity(ObjectA.GetComponent<Rigidbody>(), ObjectB.GetComponent<Rigidbody>()); //should store these GetComponents, they run too often.
            }

        }
    }

}

Hopefully, someone can make some sense of this. I'm really stuck here! Please note that I DO NOT intend to reuse 90% of this code, as it's not all mine. Once I can understand what's going on, I'll re-write for my project.

Edited April 16, 2015 by Corbald

[*Aqua*]

My comments start with // ***.

//r=norm(rVect);
//rUnitVect=rVect/r;
//v=norm(vVect);

// *** about Unity Vector3: http://docs.unity3d.com/ScriptReference/Vector3.html
// *** looks right
float r = rVect.magnitude; //matlab 'norm' is Unity Vector3.magnitude?
Vector3 rUnitVect = rVect.normalized; // rVect/r is Unity Vector3.normalized?
float v = vVect.magnitude;

//hVect=cross(rVect,vVect);
//h=norm(hVect);
//hUnitVect=hVect/h;
//ThetaUnitVect=cross(hUnitVect,rUnitVect);

Vector3 hVect = Vector3.Cross(rVect, vVect);
float h = hVect.magnitude;
Vector3 hUnitVect = hVect.normalized;
Vector3 ThetaUnitVect = Vector3.Cross(hUnitVect, rUnitVect);


//Energy=v^2/2 - muCB/r;
//sma=-muCB/(2*Energy);

float Energy = Mathf.Pow( v, 2 )/ 2 - muCB / r; //Am I doing order of operation right?
        // *** looks ok, but if you are unsure put it into brackets
        // *** example: ((v*v) / 2) - (muCB / r) - it'll calculate v*v fist, then v² / 2, then muCB / r and finally the substraction
float sma = -muCB / (2 * Energy);

//p=h^2/muCB;
//ecc=sqrt(-p/sma + 1);

float p = Mathf.Pow(h, 2) / muCB;
float ecc = Mathf.Sqrt(-p / sma + 1);

//TrueAnom=acos((p/r - 1)/(ecc));
//if(dot(rVect,vVect)<0)
//    TrueAnom=-TrueAnom;
//end

float TrueAnom = Mathf.Acos((p / r - 1) / ecc);
if (Vector3.Dot(rVect, vVect) < 0)
{
    TrueAnom = -TrueAnom;
}

//inc=acos(hUnitVect(3));

float inc = Mathf.Acos(hUnitVect[2]); //Unity vectors are from 0-2, not 1-3

//longAscNode_1=AngleZero2Pi(asin(hUnitVect(1)/sin(inc))); //is AngleZero2Pi a matlab function? what is this?
        // *** couldn't find that in the documentation, must be a custom method or from an old MATLAB version
//longAscNode_2=AngleZero2Pi(pi-asin(hUnitVect(1)/sin(inc)));
//longAscNode_3=AngleZero2Pi(acos(-hUnitVect(2)/sin(inc)));
//longAscNode_4=AngleZero2Pi(-acos(-hUnitVect(2)/sin(inc)));

//longAscNodeSet1=round(1000*[longAscNode_1,longAscNode_2])/1000; //is this a Mathf.Round as a vector2?
        // *** I think it is.
//longAscNodeSet2=round(1000*[longAscNode_3,longAscNode_4])/1000;

//[val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2); //am I finding the intersection of two vector2's? I don't know how to read this!
        // *** foo[a, b, c] is an array: 'foo' is the name of it, abc are the names of three elements inside the array.
        // *** The named elements allow direct access to array element.
        // *** Therefore [val,ia,ib] is an unnamed array with three elements, named val, ia and ib.
//longAscNode=longAscNodeSet1(ia); //I don't know how to read this!
        // *** I, too, have no idea what happens here.

//Theta_1=AngleZero2Pi(asin(rUnitVect(3)/sin(inc))); //Quesion as above
//Theta_2=AngleZero2Pi(pi-asin(rUnitVect(3)/sin(inc)));
//Theta_3=AngleZero2Pi(acos(ThetaUnitVect(3)/sin(inc)));
//Theta_4=AngleZero2Pi(-acos(ThetaUnitVect(3)/sin(inc)));

//ThetaSet1=round(1000*[Theta_1,Theta_2])/1000; //Quesion as above
//ThetaSet2=round(1000*[Theta_3,Theta_4])/1000;

//[val,ia,ib]=intersect(ThetaSet1,ThetaSet2); //Quesion as above
//Theta=ThetaSet1(ia); //Quesion as above

//ArgPeri=Theta-TrueAnom; //If I could figure out the above, this might make sense!


//some stuff will go here to draw the elipse/orbits. I'll tackle that later, after I make sense of the rest.

void FixedUpdate()
{
    //Get every object 
    GameObject[] Objects = GameObject.FindGameObjectsWithTag("Star"); //This runs too often. Extract to Start() later to optimize.
            // *** Am I right that the amount of 'Star' objects don't change while playing? In that case it might be a good idea to only
            // *** search for them once and put the result into an array.

    //the gravity between each couple of object is calculated
    foreach (GameObject ObjectA in Objects)
    {
        foreach (GameObject ObjectB in Objects)
                // *** Two nested loops are bad if there are a lot of objects because they need a runtime of (number of objects)^2.
                // *** You can do this when you know the number of objects stays small.
        {
            //Objects must not self interact 
            if (ObjectA == ObjectB)
                continue;

            ApplyGravity(ObjectA.GetComponent<Rigidbody>(), ObjectB.GetComponent<Rigidbody>()); //should store these GetComponents, they run too often.
                    // *** Same as above. If the objects don't get replaced during play, store them into a property or something like that to avoid using the methods.
        }

    }
}

If you want to understand what you are calculating then read a book about the prediction calculation of an orbit. Whoever wrote the MATLAB code did the same thing. A lot of that stuff is weird for me, too.

Edited April 17, 2015 by *Aqua*

[Corbald]

Hey Aqua, thanks for the response. Looks like I need to know how AngleZero2Pi works. Hopefully someone that understands the maths can help me figure that out. Unfortunately, I can't learn from math books, as previous experience has taught me. My mind can't stay focused to that degree. I HAVE, however, primed myself with a bit of understanding from the helpful folks who air lessons on Youtube. That said, I'm not going to understand it until I put it in practice. Thus this endeavor.

So, the questions remaining are:

How do these work?

//longAscNode_1=AngleZero2Pi(asin(hUnitVect(1)/sin(inc))); //is AngleZero2Pi a matlab function? what is this?
        // *** couldn't find that in the documentation, must be a custom method or from an old MATLAB version
//longAscNode_2=AngleZero2Pi(pi-asin(hUnitVect(1)/sin(inc)));
//longAscNode_3=AngleZero2Pi(acos(-hUnitVect(2)/sin(inc)));
//longAscNode_4=AngleZero2Pi(-acos(-hUnitVect(2)/sin(inc)));

And:

//[val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2); //am I finding the intersection of two vector2's? I don't know how to read this!
        // *** foo[a, b, c] is an array: 'foo' is the name of it, abc are the names of three elements inside the array.
        // *** The named elements allow direct access to array element.
        // *** Therefore [val,ia,ib] is an unnamed array with three elements, named val, ia and ib.
//longAscNode=longAscNodeSet1(ia); //I don't know how to read this!
        // *** I, too, have no idea what happens here.

if longAscNodeSet1 and longAscNodeSet2 are vector2's, then [val,ia,ib] should be a vector2, but it's got three elements. what the heck is 'val'? If it IS a vector 2, then longAscNodeSet1(ia) should be something like longAscNodeSet1.x, riiiight?

*Later*, I found this for [val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2):

Define two vectors with values in common.

A = [7 1 7 7 4]; B = [7 0 4 4 0];
Find the values common to both A and B, as well as the index vectors ia and ib, such that C = A(ia) and C = B(ib).

[C,ia,ib] = intersect(A,
C =

     4     7


ia =

     5
     1


ib =

     3
     1

Worst case: I write a function which does this, myself. Does unity do anything like this, built in? I'll run with the assumption that it does not (since I can't find anything in their docs) and code my own, but if anyone knows a shortcut, lemme know!

That just leaves AngleZero2Pi! Anyone have insight?

[Fel]

Hey Aqua, thanks for the response. Looks like I need to know how AngleZero2Pi works. Hopefully someone that understands the maths can help me figure that out. Unfortunately, I can't learn from math books, as previous experience has taught me. My mind can't stay focused to that degree. I HAVE, however, primed myself with a bit of understanding from the helpful folks who air lessons on Youtube. That said, I'm not going to understand it until I put it in practice. Thus this endeavor.

So, the questions remaining are:

How do these work?

//longAscNode_1=AngleZero2Pi(asin(hUnitVect(1)/sin(inc))); //is AngleZero2Pi a matlab function? what is this?
        // *** couldn't find that in the documentation, must be a custom method or from an old MATLAB version
//longAscNode_2=AngleZero2Pi(pi-asin(hUnitVect(1)/sin(inc)));
//longAscNode_3=AngleZero2Pi(acos(-hUnitVect(2)/sin(inc)));
//longAscNode_4=AngleZero2Pi(-acos(-hUnitVect(2)/sin(inc)));

And:

//[val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2); //am I finding the intersection of two vector2's? I don't know how to read this!
        // *** foo[a, b, c] is an array: 'foo' is the name of it, abc are the names of three elements inside the array.
        // *** The named elements allow direct access to array element.
        // *** Therefore [val,ia,ib] is an unnamed array with three elements, named val, ia and ib.
//longAscNode=longAscNodeSet1(ia); //I don't know how to read this!
        // *** I, too, have no idea what happens here.

if longAscNodeSet1 and longAscNodeSet2 are vector2's, then [val,ia,ib] should be a vector2, but it's got three elements. what the heck is 'val'? If it IS a vector 2, then longAscNodeSet1(ia) should be something like longAscNodeSet1.x, riiiight?

*Later*, I found this for [val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2):

Define two vectors with values in common.

A = [7 1 7 7 4]; B = [7 0 4 4 0];
Find the values common to both A and B, as well as the index vectors ia and ib, such that C = A(ia) and C = B(ib).

[C,ia,ib] = intersect(A,
C =

     4     7


ia =

     5
     1


ib =

     3
     1

Worst case: I write a function which does this, myself. Does unity do anything like this, built in? I'll run with the assumption that it does not (since I can't find anything in their docs) and code my own, but if anyone knows a shortcut, lemme know!

That just leaves AngleZero2Pi! Anyone have insight?

http://www.google.com/search?q=%2BAngleZero2Pi&spell&hl=en&sa=X&&as_q&nfpr=1

Custom function, my guess? some kind of wrap.

Angle = 540;

Equivalent Angle = 540 - 360 = 180

http://www.mathworks.com/help/map/ref/wraptopi.html

longAscNodeSet1=round(1000*[longAscNode_1,longAscNode_2])/1000;
longAscNodeSet2=round(1000*[longAscNode_3,longAscNode_4])/1000;

[val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2);
longAscNode=longAscNodeSet1(ia);

this just says, stupidly (the performance gain for using SSE code is near zero) [actually, I would argue that it is worse off because it is bad SSE code]

if (longAscNode_1 == longAscNode_3 || longAscNode_1 == longAscNode_4) longAscNode = longAscNode_1;
if (longAscNode_2 == longAscNode_3 || longAscNode_2 == longAscNode_4) longAscNode = longAscNode_2;

Of course, "if both" then you get a bug where longAscNode = [longAscNode_1, longAscNode_2]

Edited April 17, 2015 by Fel

[Corbald]

Fel, Aqua, you two are CHAMPS! Chuck Norris aspires to be like you! This should be enough for me to get something working!

Fel, I should have realized that the longAscNode gobbledygook was if's. Makes sense when you run numbers through it! DOH! The Lambda wrapping function would have kept me guessing for weeks! I kept getting lost in the interior of the function's argument calls! Thanks for clearing that up, especially.