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)

Question

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 prograde = true;
    const lowPath = true;
    const vals = await solver(
      mu,
      r0.invertY().toArray(),
      r1.invertY().toArray(),
      deltat,
      0,
      prograde,
      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.add(v1.invertY().multiply(2000000)),
      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
    const secondFocus = center.add(direction.multiply(c));

    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
  public add(v: Vector): Vector {
    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.

Answer 1

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);
    const p = c1.add(c2.subtract(c1).normalize().multiply(a));

    const x1 = p.x + (h * (c2.y - c1.y)) / d;
    const y1 = p.y - (h * (c2.x - c1.x)) / d;

    const x2 = p.x - (h * (c2.y - c1.y)) / d;
    const y2 = p.y + (h * (c2.x - c1.x)) / d;

    return [new Vector(x1, y1), new Vector(x2, y2)];
  }

resulting in these ellipses, including the one I need! It also matches up perfectly with the initial velocity vector, further confirming it's the result I need.