# How do i correctly calculate the ellipse that illustrates the orbit of an object in space based on initial state values (r0, v0, deltat, and mu)

I'm working on a game that has a primary UI showing the solar system in a top-down (earth's north pole being roughly 'up'). And I have a need to show a semi-realistic path of travel for an object between planets. I have been using Lambert's problem to determine the initial velocity vector for the object, and the results I'm getting from that seem reasonable for my game.

I've been trying to solve this for a couple weeks now, and have even tried using ChatGPT-4 to help with the math that I don't understand. I have come close a couple of times, but never seem to be able to land at the appropriate solution.

My calculations for a (semi-major-axis of the ellipse) sometimes come up negative, which I know is incorrect.

One cause that I considered, is my game (like most) use a left-handed coordinate system, but most of the math expects right-handed. (in my game, +y is down, in most of the math +y is up). I've dealt with this by inverting my y value any time I do calculations on the vectors, and making sure the result is inverted any time I'm displaying to the screen. I have also tried inverting my entire coordinate system, but the results were similar.

below is the TS implementation of the logic generating my current ellipse.

import { Constants } from './constants';
import Vector from './vector';
import * as math from 'mathjs';
const solver = require('lambert-orbit');

export default class LambertSolver {
public static async generatePath(
r0: Vector, // Initial position vector (km)
r1: Vector, // Final position vector (km)
deltat: number, // Time of flight (s)
solarMass: number, // Solar mass (kg)
pathResolution: number = 200
): Promise<Vector[]> {
const mu = solarMass * Constants.GRAVITY_KM;
const lowPath = true;
const vals = await solver(
mu,
r0.invertY().toArray(),
r1.invertY().toArray(),
deltat,
0,
lowPath
);

const v1 = Vector.fromArray(vals[0]);

const [a] = this.orbitalParameters(r0.invertY(), v1, mu);

// const v1Norm = v1.magnitude();
// const r1Norm = r1.magnitude();

// const a = -mu / (v1Norm * v1Norm - (2 * mu) / r1Norm);

const secondFoci = this.findSecondFoci(r0.invertY(), v1, a, mu);
// const [secondFoci2, secondFoci3] = this.findSecondFociAlternative(
//   r0.invertY(),
//   r1.invertY(),
//   a
// );

const path = this.ellipsePath(
new Vector(0, 0),
secondFoci,
a / 0.9,
pathResolution
);

return [
...path.map((x) => x.invertY()),
// ...this.ellipsePath(new Vector(0, 0), secondFoci2, a, pathResolution).map(
//   (x) => x.invertY()
// ),
// ...this.ellipsePath(new Vector(0, 0), secondFoci3, a, pathResolution).map(
//   (x) => x.invertY()
// ),
r0,
r1,
secondFoci.invertY(),
// secondFoci2,
// secondFoci3,
new Vector(0, 0),
r0,
];
}

private static orbitalParameters(
r1: Vector,
v1: Vector,
mu: number
): [number, number, number, number, number, number] {
const R = r1.magnitude();
const V = v1.magnitude();
const a = 1 / (2 / R - (V * V) / mu);

const Hv = r1.cross(v1);
const H = Hv.magnitude();

const p = (H * H) / mu;

const q = r1.dot(v1);

const e = Math.sqrt(1 - p / a);

const i = Math.acos(Hv.z / H);

let omega = 0;
if (i != 0) {
omega = Math.atan2(Hv.x, -Hv.y);
}
const Cw = (r1.x * Math.cos(omega) + r1.y * Math.sin(omega)) / R;

let Sw = 0;
if (i == 0 || i == Math.PI) {
Sw = (r1.y * Math.cos(omega) - r1.x * Math.sin(omega)) / R;
} else {
Sw = r1.z / (R * Math.sin(i));
}
const TAx = (H * H) / (R * mu) - 1;
const TAy = (H * q) / (R * mu);
const v = Math.atan2(TAy, TAx);

var w = Math.atan2(Sw, Cw) - v;

return [a, e, i, omega, w, v];
}

private static findSecondFoci(
r1: Vector,
v1: Vector,
a: number,
mu: number
): Vector {
// Compute the specific angular momentum vector
const h = r1.cross(v1);

// Compute the eccentricity vector
const epsilon = v1.magnitude() ** 2 / 2 - mu / r1.magnitude();
const e_vec = r1
.multiply(epsilon)
.subtract(v1.multiply(r1.dot(v1)))
.multiply(1 / mu);

// Compute the eccentricity (magnitude of the eccentricity vector)
const e = e_vec.magnitude();

// Compute the distance from the center of the ellipse to each focus
const c = a * e;

// Compute the center of the ellipse
const firstFocus = r1;
const center = firstFocus.multiply(1 - e);

// Find the direction vector from the first focus to the second focus
const direction = e_vec.normalize();

// Compute the position of the second focus

return secondFocus;
}

private static ellipsePath(
focalPoint1: Vector,
focalPoint2: Vector,
a: number,
numberOfPoints: number = 50
): Vector[] {
// Calculate the distance between the two focal points
const focalDistance = focalPoint1.distanceTo(focalPoint2);

// Calculate the semi-minor axis (b)
const c = focalDistance / 2;
const b = Math.sqrt(a * a - c * c);

// Calculate the center of the ellipse
const centerX = (focalPoint1.x + focalPoint2.x) / 2;
const centerY = (focalPoint1.y + focalPoint2.y) / 2;
const center = new Vector(centerX, centerY);

// Calculate the angle of the line connecting the focal points with the positive x-axis
const angle = Math.atan2(
focalPoint2.y - focalPoint1.y,
focalPoint2.x - focalPoint1.x
);

// Generate ellipse points in polar coordinates and convert them to Cartesian coordinates
const path: Vector[] = [];
for (let i = 0; i < numberOfPoints; i++) {
const theta = (2 * Math.PI * i) / numberOfPoints;

// Ellipse equation in polar coordinates: r(theta) = (a * b) / sqrt((b * cos(theta))^2 + (a * sin(theta))^2)
const r =
(a * b) /
Math.sqrt(
Math.pow(b * Math.cos(theta), 2) + Math.pow(a * Math.sin(theta), 2)
);

// Convert polar coordinates to Cartesian coordinates
const x = r * Math.cos(theta);
const y = r * Math.sin(theta);

// Rotate the point by the angle and translate it to the center
const rotatedX = x * Math.cos(angle) - y * Math.sin(angle);
const rotatedY = x * Math.sin(angle) + y * Math.cos(angle);

const translatedX = rotatedX + centerX;
const translatedY = rotatedY + centerY;

path.push(new Vector(translatedX, translatedY));
}

return path;
}
}


After having several version of Lambert problem solvers in my own code (AI assisted), I settled on using the lambert-orbit package, which seems to do exactly what I need, and adds in some checks and additional logic I'll need later.

My vector class implementation is below:

export default class Vector {
public x: number;
public y: number;
public z: number;

constructor(x: number, y: number, z: number = 0) {
this.x = x;
this.y = y;
this.z = z;
}

// Add another vector to this one
return new Vector(this.x + v.x, this.y + v.y, this.z + v.z);
}

// Subtract another vector from this one
public subtract(v: Vector): Vector {
return new Vector(this.x - v.x, this.y - v.y, this.z - v.z);
}

// Multiply the vector by a scalar
public multiply(scalar: number): Vector {
return new Vector(this.x * scalar, this.y * scalar, this.z * scalar);
}

// Divide the vector by a scalar
public divide(scalar: number): Vector {
return new Vector(this.x / scalar, this.y / scalar, this.z / scalar);
}

// Calculate the magnitude of the vector
public magnitude(): number {
return Math.sqrt(this.x * this.x + this.y * this.y + this.z * this.z);
}

// Calculate the dot product of two vectors
public dot(v: Vector): number {
return this.x * v.x + this.y * v.y + this.z * v.z;
}

public distanceTo(other: Vector): number {
const dx = this.x - other.x;
const dy = this.y - other.y;
const dz = this.z - other.z;
return Math.sqrt(dx * dx + dy * dy + dz * dz);
}

// Copy the components of another vector into this one
public copy(v: Vector): Vector {
this.x = v.x;
this.y = v.y;
this.z = v.z;
return this;
}

public scale(scalar: number): Vector {
return new Vector(this.x * scalar, this.y * scalar, this.z * scalar);
}

public descale(scalar: number): Vector {
return new Vector(this.x / scalar, this.y / scalar, this.z / scalar);
}

public unit(): Vector {
const magnitude = this.magnitude();
if (magnitude === 0) {
throw new Error('Cannot compute the unit vector of a zero vector.');
}
return new Vector(
this.x / magnitude,
this.y / magnitude,
this.z / magnitude
);
}

public absolute(): Vector {
return new Vector(Math.abs(this.x), Math.abs(this.y), Math.abs(this.z));
}

// Create a new vector from this one
public clone(): Vector {
return new Vector(this.x, this.y, this.z);
}

public toArray(): number[] {
return [this.x, this.y, this.z];
}

public static fromArray(array: number[]): Vector {
return new Vector(array[0], array[1], array[2]);
}

public angleTo(other: Vector): number {
const dx = other.x - this.x;
const dy = other.y - this.y;
const dz = other.z - this.z;
const d = Math.sqrt(dx * dx + dy * dy + dz * dz);
return Math.acos(this.dot(other) / (this.magnitude() * d));
}

public invert(): Vector {
return new Vector(-this.x, -this.y, -this.z);
}

public invertY(): Vector {
return new Vector(this.x, -this.y, this.z);
}

public invertX(): Vector {
return new Vector(-this.x, this.y, this.z);
}

public invertZ(): Vector {
return new Vector(this.x, this.y, -this.z);
}

public angle(): number {
return Math.atan2(this.y, this.x);
}

public normalize(): Vector {
const norm = this.magnitude();
return new Vector(this.x / norm, this.y / norm, this.z / norm);
}

public cross(other: Vector): Vector {
const x = this.y * other.z - this.z * other.y;
const y = this.z * other.x - this.x * other.z;
const z = this.x * other.y - this.y * other.x;
return new Vector(x, y, z);
}
}


For an example, if I pass in the following to the generatePath method: r0: {x: 126619585.3943946, y: -79668083.39934838} r1: {x: 78468331.12907715, y: -214058085.17912537} deltat: 7776000 solarMass: 1.989e+30

I get a velocity vector of: {x: 4.544347796414998, y: 28.198168671041273}

My orbitalParameters method returns an a value of 138428143.17335743

And my findSecondFoci method returns {x: -82494511.2804895, y: -126466117.64851947}

Which results in the following screenshot. Sorry, it's a little bit of a mess, but I just started adding points to the drawn path for debugging. the screenshot has descriptions of what each point are visualized in the game space.

As you can see, the ellipse is not only not even close to intercepting both r0 and r1 locations, it seems to be in the wrong direction.

Interestingly after seeing this, I've taken another screenshot showing the X values of the path inverted instead of Y, and the tip of the ellipse is remarkably close to the target destination. Doing some testing now to see if this is just a coincidence. If this is not a coincidence, I just need to figure out why my semi-minor-axis is so large.

• As a note - the second image seems to have been a coincidence. It doesn't look anything like this beyond that point in the orbit Apr 25, 2023 at 15:39
• Never use chatgpt if you need a correct answer. Apr 25, 2023 at 18:22
• It's helped me significantly in many cases, even helped me arrive at a working lambert solver, but it did take a lot of time and effort. gpt-3 is nowhere near as good as gpt-4. In either case, I just don't have the fundamental algebra/trig knowledge needed to come to my own solution. @OrganicMarble do you have any suggestions? Apr 25, 2023 at 18:48

Finally found the solution I was looking for!

I came across this graphic showing the visualization of the ellipse.

From That I generated some code to find me the two possible locations for the second focus point, by generating a circle that represents 2a - r and 2a - r0 around each of their respective points. Where they intersect, are my possible focal points.

Here's the code:

private static circleCircleIntersection(
c1: Vector,
r1: number,
c2: Vector,
r2: number
): [Vector, Vector] | null {
const d = c2.subtract(c1).magnitude();
if (d > r1 + r2 || d < Math.abs(r1 - r2)) {
return null;
}

const a = (r1 * r1 - r2 * r2 + d * d) / (2 * d);
const h = Math.sqrt(r1 * r1 - a * a);