How to navigate a moving 2d spaceship under newtonian mechanics to a point in space?

Question

I am creating a real-time game where ships move through space using basic laws of motion. The ships move by rotating (instantly for our purposes) to an angle and then firing their engine to accelerate.

I have coded up a basic ECS system wherein every tick velocity is added to position and acceleration is added to velocity.

So the problem is at any given moment how can a ship determine the rotation to fire its engine (assume constant acceleration) to intercept an arbitrary point elsewhere in the world. Ideally in the shortest time possible, however I'd settle for ever.

I have the following variables:\$v\$ the current velocity vector of the ship at \$t = 0\$ , \$P_0\$ the current postition of the ship, \$P_2\$ the desired intercept point, and \$||a||\$ the maximum magnitude of the acceleration vector.

I need to find some rotation \$\theta\$ to fire my engine in to follow a curve to the given intercept point \$P_2\$.

I have tried solving for the quadratic bezier curve control points for such a path and then deriving its acceleration. Like so:

\$P'(t) = 2(1-t)(P_1 - P_0)+2t(P_2 - P_1)\$

\$P'(0) = v = 2P_1 - 2P_0\$

\$P_1 = \frac{v}{2} + P_0\$

To derive \$P_1\$ and then:

\$P''(t) = 2(P_2 - 2P_1 + P_0)\$

plugging in the three points to get an acceleration vector, however if I understand correctly bezier curves are from \$t=0\$ to \$t=1\$? So this is a massive acceleration vector that would get my ship to the point in one game tick. I have tried scaling this vector to the maximum magnitude of the engine's acceleration, but that doesn't quite work. I feel like my approach is close but I can't quite figure it out.

Cont.

I have thought about this more and returned to this solution. I believe I need to solve for time \$T\$ to reach the point by calculating the radius of the circle in the aforementioned answer. If I do that then I can calculate the second control point as \$P_0 + Tv\$.

It seems to me that the minimum time necessary to reach a point is when the destination point \$P_2\$ is on the circumference of the circle described by it's center \$P_0 + Tv\$ and radius \$\frac{T^2}{2}a\$ where \$a\$ is the maximum acceleration.

Previously I was conflating the \$t\$ in the bezier curve polynomial (between 0 and 1) with some time \$T\$ it will take to complete the maneuver which is unbounded.

Answer 1

I found a working solution for myself. Here's my code in Typescript (using bitecs) with comments:

// eid here is the entity id in my ECS
const position = getPos(eid);
const velocity = getVelocity(eid);
const destination = getDestination(eid);
const mass = Mass.value[eid];
const maxForce = Engine.maxForce[eid];

// Bezier curve control points
// p0 and p2 are trivial: they're just the starting and ending location
// p1 is half the distance the ship would have travelled if no further
// acceleration. The core issue is we must know the entirety of the
// travel time in terms of the unit our velocity is in.
// In this case frames/ticks.
const p0 = position;
// this function is explained below. It returns the time used for the 
// entire maneuver
const time = findTimeToPoint(p0, destination, maxForce / mass, velocity);
// Here we add the distance to p0 that is travelled in time / 2 frames
const p1 = vec2Add(p0, vec2Mult(velocity, time / 2));
const p2 = destination;

// The second derivative of a quadratic Bezier curve is
// P''(t) = 2(P_2 - 2P_1 + P_0)
// which makes it a constant value with respect to time. Since we have
// all three points we can just solve for it.
const a = vec2Mult(vec2Add(vec2Sub(p2, vec2Mult(p1, 2)), p0), 2);

// The resulting vector (a) of the acceleration has a large magnitude
// since it's in terms of unitary time from 0 to 1. However since we only
// need the rotation and we already know the magnitude we can just
// solve it.
let rotation = Math.atan(a.y / a.x);
// deal w/ principle values of atan
if (a.x < 0) {
  rotation += Math.PI;
}

// Then we update the rotation and my physics systems take it from there
Engine.percent[eid] = 1;
Rotation.radians[eid] = rotation;
removeComponent(world, MoveTo, eid);

Here is my function to find the time of the entire maneuver. It's not very efficient but works for my use case. You could potentially solve for T using the proper equations:

$$ \displaylines{ S = \|\vec{P_3} - \vec{v}T\| \\ S = \sqrt{(P_{3x} - v_xT)^2 + (P_{3y} - v_yT)^2} } $$

and

$$ R = \frac{T^2}{2}a $$

but I couldn't figure it out. This sort of computational iterative approach is good enough for now.

const findTimeToPoint = (p0: Vec2, p3: Vec2, a: number, v: Vec2) => {
  // Coordinate transform in terms of origin at p1
  // I found it easier to think about p3 moving _away_ from 
  // the center of the circle than the circle moving but I think
  // it's the same
  const p3_ = vec2Sub(p0, p3);
  const v_ = vec2Mult(v, -1);

  let T = 0;
  // The distance between the potential position without acceleration
  // at time T and the edge of a circle described by R
  let S = vec2Magnitude(p3_);
  let R = 0;
  // We just iterate one frame forward and solve for R and S
  // when R < S then the time given is enough for the movement
  // from p0 to p3
  while (R < S) {
    T++;
    S = vec2Magnitude(vec2Sub(p3_, vec2Mult(v_, T)));
    R = a * (Math.pow(T, 2) / 2);
  }
  return T;
};

This works to have the ship follow a curve to an intercept point. I can probably do some cool stuff with the engine percentage if I wanted to take fuel usage into account but this is good enough for now.