Calculating a collision course

[ThePuska]

You're the navigation computer of a space ship. You've been fighting a superior enemy for some time, and now your weapon systems are damaged, your shields are failing and you're leaking atmosphere. As his last, desperate attempt at damaging the enemy, your captain orders you: "Take us on a collision course!" What do you do? What is a collision course to a computer? How do you calculate it?

Most significantly it's one of many solutions to a mathematical equation that equates two trajectories: the projectile's and the target's. But it's not just any solution, it's one that minimizes the time until collision and maximizes the crashing speed with the target.

[img]http://img155.imageshack.us/img155/1838/whycollide.jpg[/img]

So the problem is just a mathematical equation. A computer can understand that.

Our target is any polynomial vector function of time. It can be another space ship flying in a straight path, or accelerating, or changing its acceleration predictably. It can be an object orbiting another. Though an orbit can't be exactly expressed as a polynomial function, a part of its curve can be approximated as such.

The target function absorbs our own displacement, initial velocity and gravitational acceleration. Our ship's function is left only with our own acceleration; the one our ship's engines are responsible for.

[release][h2]The Math[/h2]

So let's solve the collision course. We know how the target is moving. We know what kind of thrust our own ship is capable of producing. We just need to know how to utilize this to crash into the target.

Let the target function any polynomial with vector coefficients
[img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{F}(t)%20=%20%5Csum%20_{k=0}%20^N%20{%5Coverrightarrow{a_k}%20t^k}[/img]

Our projectile's function is
[img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{G}(t)%20=%20%5Cint%20%5Cint%20%5Coverrightarrow{A}(x)dxdx[/img]
i.e. the second integral of an acceleration function. Assume that the constants that the integrals require are zero.

We'll simplify it by only considering constant acceleration. When you want to crash into something, you don't slow down. This gives us a collision function of
[img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{G}(t)%20=%20{1%20%5Cover%202}%5Coverrightarrow{A}t^2[/img]

Writing these trajectories equal gives us
[img]http://latex.codecogs.com/gif.latex?%5Csum%20_{k=0}%20^N%20{%5Coverrightarrow{a_k}%20t^k}%20=%20{1%20%5Cover%202}%5Coverrightarrow{A}t^2[/img]

Which isn't very useful on its own. We don't know the direction of our acceleration. Is it directly towards the target? Slightly off of it?

However if we assume that our ship's engine is pushed to its full capabilities for the whole duration of the acceleration, we can rewrite the right side as a product of the maximum acceleration our engine is capable of producing and a direction. The direction is unknown to us, but we know our own maximum acceleration.

[img]http://latex.codecogs.com/gif.latex?%5Csum%20_{k=0}%20^N%20{%5Coverrightarrow{a_k}%20t^k}%20=%20{1%20%5Cover%202}T%5Coverrightarrow{d}t^2[/img]

We can get rid of the direction vector by squaring both sides, as the square of a unit vector is unity. [url=http://www.facepunch.com/threads/1091733-Calculating-a-collision-course?p=30060164&viewfull=1#post30060164](Clarification)[/url]

[img]http://latex.codecogs.com/gif.latex?(%5Csum%20_{k=0}%20^N%20{%5Coverrightarrow{a_k}%20t^k})^2%20=%20{1%20%5Cover%204}T^2t^4[/img]

Now this is more useful. [i]The only unknown variable in this equation is t, so we can solve it.[/i] As a polynomial equation it's easy enough. As long as the target trajectory is quadratic or simpler, the equation can be solved with some symbolic manipulation. For more generic equations, numerical analysis is the way to go.[/release]

[img]http://img854.imageshack.us/img854/1067/stargate.jpg[/img]

[release][h2]The Programming[/h2]

This requires some vector and complex number math. I'm not going to explain the theory behind those, but the bare minimum implementations that are required can be found here:
[url=http://pastebin.com/KU9G1xAs]ComplexDef.cs[/url]
[url=http://pastebin.com/HWyhiciy]Vectors.cs[/url]

It's worth noting that only basic arithmetic is needed. I also used a function for raising complex numbers to integer powers using repeated multiplication. This should not be replaced with more generic exponentiation functions as not all exponentiation is equivalent in the complex plane.

The math works in any number of dimensions but I limited my vectors to three dimensions.

The polynomial functions are expressed simply as an array of doubles or vectors. Each double or vector corresponds to a term's coefficient, and the term's degree is directly its position in the array. The zeroth index thus corresponds to a constant.

[img]http://img121.imageshack.us/img121/7338/polyinmemory.png[/img]

Thus the function that is used to calculate a collision course with a given target trajectory has the signature:

[csharp]public static CollisionResult CalculateCollision(Vec3d[] target, double maxAcceleration);[/csharp]

Where the struct CollisionResult contains:

[csharp]public struct CollisionResult
{
public bool Collision; // Can a collision take place?
public double Time; // How long to the optimal collision?
public Vec3d Acceleration; // How do we need to accelerate to collide?
public Vec3d Position; // Where do we collide?
}[/csharp]

The flow of the function is:

[list]
[*][b]Construct a polynomial function (with real coefficients) of the target function and our maximum acceleration.[/b]

[*][b]Solve every root of the polynomial function. Each root corresponds to a possible collision.[/b]

[*][b]Pick the best collision.[/b]

[*][b]Calculate and return information about the collision.[/b]
[/list]

Let's step through the function in parts.

The first thing to do is to turn the last equation of the math part into code. We'll throw it all on the left side of the equation, leaving us a polynomial function whose coefficients are all real numbers (with vector multiplication meaning the dot product).

[img]http://latex.codecogs.com/gif.latex?(%5Csum%20_{k=0}%20^N%20{%5Coverrightarrow{a_k}%20t^k})^2%20-%20{1%20%5Cover%204}T^2t^4%20=%200[/img]

With [img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{a_0}%20..%20%5Coverrightarrow{a_N}[/img] corresponding to the vectors in [i]target[/i] and T being equal to [i]maxAcceleration[/i], we have:

[csharp]// Construct the polynomial function

// This means squaring the target function and subtracting from it the square of the acceleration function

// The maximum degree is either (target function's highest degree)^2 or (our acceleration's degree)^2
double[] poly = new double[Math.Max(2 * (target.Length - 1), 4) + 1];
{
int i, j;

for (i = 0; i < target.Length; i++)
{
for (j = 0; j < target.Length; j++)
{
// The full terms, when written out with their degrees, are
// (target[i] * t^i) * (target[j] * t^j) = (target[i] * target[j] * t^(i + j))
poly[i + j] += Vec3d.DotP(target[i], target[j]);
}
}

poly[4] -= maxAcceleration * maxAcceleration / 4;
}
[/csharp]

After that we need to solve the polynomial function. Any method will work.

[csharp]// Solve the polynomial function

// We'll use the Durand-Kerner method.
// We don't actually need the polynomial function after this so we'll just change it a bit.
// The D-K method requires it to be monic, so we'll divide every term by biggest term's coefficient.

for (int i = 0; i < poly.Length; i++)
{
poly[i] /= poly[poly.Length - 1];
}

Complex[] polySolutions = PolyRoots(poly);[/csharp]

The function PolyRoots has the signature
[csharp]public static Complex[] PolyRoots(double[] poly)[/csharp]

and it can be found here: [url]http://pastebin.com/iKXUS2FB[/url]

Meaningless roots are culled and the best one is picked

[csharp]// The solutions are in no particular order and contain every complex root of the polynomial.
// We don't know what imaginary time is and we're not interested in negative time. Cull useless roots.

const double imTreshold = 0.01; // This is an arbitrary constant
int rootCount = polySolutions.Length;
for (int i = polySolutions.Length - 1; i >= 0; i--)
{
if ((polySolutions[i].Re < 0) || (Math.Abs(polySolutions[i].Im) > imTreshold))
{
// Delete this root
polySolutions[i] = polySolutions[--rootCount];
}
}

// If that already exhausted the root pool, return nothing useful

if (rootCount == 0)
{
return new CollisionResult(false, 0, new Vec3d(0, 0, 0), new Vec3d(0, 0, 0));
}

// Otherwise search for the best root

int bestRoot = 0;
double time = polySolutions[0].Re;
for (int i = 1; i < rootCount; i++)
{
if (polySolutions[i].Re < time)
{
bestRoot = i;
time = polySolutions[i].Re;
}
}[/csharp]

After which it's simple to get more information out of the collision.

[csharp]// Calculate the acceleration & position with this time

Vec3d position = new Vec3d(0, 0, 0);
for (int i = 0; i < target.Length; i++)
{
position += Math.Pow(time, i) * target[i];
}
Vec3d acceleration = (2 / (time * time)) * position;

return new CollisionResult(true, time, acceleration, position);[/csharp]

And that's it![/release]

[release][b]The whole code:[/b]

[url=http://pastebin.com/BNpEnrzK]Solver.cs[/url]
[url=http://pastebin.com/KU9G1xAs]ComplexDef.cs[/url]
[url=http://pastebin.com/HWyhiciy]Vectors.cs[/url]
[/release]

[release][h2]How is it used?[/h2]

The input vector polynomial needs to describe the motion of the target as a function of time. It could be based on observation, or other knowledge of the target's motion.

For example, on Earth, a thrown object has the trajectory
[img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{F}(t)%20=%20%5Coverrightarrow{x_0}%20+%20%5Coverrightarrow{v_0}t%20+%20{1%5Cover2}%5Coverrightarrow{g}t^2[/img]

Where [img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{x_0}[/img] is the location of the thrower, [img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{v_0}[/img] is the object's initial velocity and [img]http://latex.codecogs.com/gif.latex?%5Coverrightarrow{g}[/img] is the gravitational acceleration.

The corresponding vector array that you would pass to the function is [b][x0, v0, 0.5*g][/b]

As a result you will get a prediction of when the collision will happen, where it will happen and what kind of constant acceleration your engine must produce in order for it to happen.

What it can't take into account is things that do not depend directly on time, but on something else, like the object's position in space. As an example, a changing gravitational acceleration, in the case of very long distances, will have to be dealt with somehow else.
[/release]

[Anonim]

I can tell this is going to be good as a reference thread later.

[Irtimid]

I think the derivative of acceleration is called "jerk."

[Remscar]

I wouldnt crash my ship into their mothership. I would call the asgard to fuck them up.

[benjojo]

Oddly I just learnt this in physics.

[ThePuska]

[QUOTE=ThePuska;30047810]We can get rid of the direction vector by squaring both sides, as the square of a unit vector is unity.

[img]http://latex.codecogs.com/gif.latex?(%5Csum%20_{k=0}%20^N%20{%5Coverrightarrow{a_k}%20t^k})^2%20=%20{1%20%5Cover%204}T^2t^4[/img][/QUOTE]

This might need some additional explanation, as it isn't always a valid method.

It relies on the implication that
[img]http://latex.codecogs.com/gif.latex?f(x_0)%20=%20f(x_1)[/img]
means
[img]http://latex.codecogs.com/gif.latex?x_0%20=%20x_1[/img]

Which is only true if f(x) is injective. The function we used, squaring, is not injective. Clearly, if you have
[img]http://latex.codecogs.com/gif.latex?x%20=%205[/img]
and you square it
[img]http://latex.codecogs.com/gif.latex?x^2%20=%2025[/img]
you can get two solutions for x, only one of which is correct.
[img]http://latex.codecogs.com/gif.latex?x%20=%20%5Cpm%205[/img]

This can be avoided by only considering a part of its domain at a time. We're only interested in positive time in the actual implementation. The square function is strictly increasing with positive input, which makes it injective, and the implication above holds.

[Zwolf11]

[img]http://dl.dropbox.com/u/7958512/teehee.png[/img]

Teehee, it looks like a face

[icantread49]

You use tons of math yet you depict speed as a vector :eng99:

[ThePuska]

[QUOTE=icantread49;30070762]You use tons of math yet you depict speed as a vector :eng99:[/QUOTE]

I'd be glad if that was the only mistake I made