I am using Unity engine and C#, and am trying to work out a way to plot future 2D orbital trajectories. I currently have a two-body setup. The planet has a mass of 1000000, the satellite a mass of 1.
In the image above you can see a pink line indicating the orbital path (PAST, NOT FUTURE) created by a trail renderer, which indicates the math is wrong since the eccentricity is not anywhere between 0 and 1.
I have also noticed that the AP and PE are incorrect since the PE is showing up as -AP. I know my math is wrong somewhere but have no idea where.
If someone much smarter than I is able to figure it out/fix this, it would be great if they could also point me in the direction of how I would go about using this to plot a future trajectory in real time.
void JustKeplerThings(Vector3 curPos, Vector3 shipCurVel, GameObject orbiting) {
//Gravitational Constant
float G = 6.67408f * Mathf.Pow(10, -11);
//Mass of ship
float m1 = thisShip.mass;
//Mass of the larger body
float m2 = -orbiting.GetComponent<PointEffector2D>().forceMagnitude;
float M = m2;
//Standard Gravitational Parameter
float µ = G * (m1 + m2);
//float µ = G * M;
//Relative Position Vector
Vector3 r = FindRelativePosition(transform, orbitee.transform.position);
//Relative Velocity Vector
Vector3 v = shipCurVel - new Vector3(orbitee.GetComponent<Rigidbody2D>().velocity.x, orbitee.GetComponent<Rigidbody2D>().velocity.y, 0);
//Specific Angular Momentum
Vector3 h = Vector3.Cross(r, v);
//Eccentricity Vector
Vector3 evec = (Vector3.Cross(v, h) / µ) - (r / Vector3.Magnitude(r));
//Eccentricity
float e = Vector3.Magnitude(evec);
//Mechanical Energy
float E = (Mathf.Pow(Vector3.Magnitude(v), 2) / 2) - (µ / Vector3.Magnitude(r));
//Semi-Major Axis
float a = µ / (2 * E);
if (Mathf.Abs(e - 1.0f) <= Mathf.Epsilon)
a = Mathf.Infinity;
//Semi-Latus Rectum Check
float p = a * (1 - Mathf.Pow(e, 2));
if (Mathf.Abs(e - 1.0f) <= Mathf.Epsilon)
p = Mathf.Pow(Vector3.Magnitude(h), 2) / µ;
//Inclination (2D)
float i = 0;
//Longitude of Ascending Node (2D)
float Ω = 0;
//Argument of Periapsis
float ω = Mathf.Atan2(evec.y, evec.x);
//ω = (2 * Mathf.PI) - ω;
//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;
//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) / µ);
}
-
public static Vector3 FindRelativePosition(Transform origin, Vector3 position) {
Vector3 distance = position - origin.position;
Vector3 relativePosition = Vector3.zero;
relativePosition.x = Vector3.Dot(distance, origin.right.normalized);
relativePosition.y = Vector3.Dot(distance, origin.up.normalized);
relativePosition.z = Vector3.Dot(distance, origin.forward.normalized);
return relativePosition;
}