Calculating an intercept course to a target with constant direction and velocity (in a 2-dimensional plane)


Warning: WP_Syntax::substituteToken(): Argument #1 ($match) must be passed by reference, value given in /home/www/jaran.de/goodbits/wp-content/plugins/wp-syntax/wp-syntax.php on line 380

For some game idea I worked on a bit some time ago I needed a bit of maths. Luckily all that was needed was taught at school and with a bit of thinking (and trying) it all worked out:
diagram of a moving object to be intercepted by a projectile

Let’s consider an object at point A moving with the constant speed vector \vec{v}. At point B we got the ability to fire a projectile with speed s. Where do we have to aim at to hit the object with our projectile? In other words: what is the location of point I or what would be the vector \vec{x}?

Now let us derive the necessary maths and put it into code… (requires understanding of vectors and quadratic equations)

If t denotes the time (such that at t=0 we got our object at point A and our “gun” at point B then the point of interception will be I = \vec{a}+t\cdot\vec{v}. We also have I = \vec{b}+t\cdot\vec{x} since we’re trying to find the point at which object and projectile are at the same place (I) at the same time (t). If we can calculate t we can also calculate I.

A first attempt would be to set equal our two (equal) locations of I that is \vec{a}+t\cdot\vec{v} = \vec{b}+t\cdot\vec{x}. However, we don’t know vector \vec{x} but only its length (i.e. the speed of our projectile) |\vec{x}| = s. Using that last information you can set up an equation system and try to solve it. It’s possible, but very nasty, so let’s take a different approach:

Since the projectile needs time t to get from B to I we can divide this distance by the known speed s to just derive that t. The distance from B to I is |\vec{a}+t\cdot\vec{v} - \vec{b}| and thus we get t=\frac{|\vec{a}+t\cdot\vec{v} - \vec{b}|}{s}. Now we got the problem that our t is on both sides of the equation which unfortunately cannot be solved by “simple” means like multiplying and dividing both sides of the equation with whatever. It appears we’ve got to get rid of the vectors first but before doing so let’s clean up things a bit:

The numerator contains the term \vec{a} - \vec{b} which are two known figures in our problem. Hence we will just use \vec{o} (o like offset between a and b) instead, that is \vec{o} = \vec{a} - \vec{b} which gives us the simpler equation t=\frac{|\vec{o}+t\cdot\vec{v}|}{s}. To get rid of the fraction we also multiply both sides with s and get s \cdot t = |\vec{o} + t \cdot \vec{v}|. Now that’s better! Just let’s get rid of the vectors and their vector length operator in one fell swoop.

Since the length of a vector is the square root of the sum of the squares of its components we get s \cdot t = \sqrt{(o_x + t \cdot v_x)^2+(o_y + t \cdot v_y)^2}, where the subscripts denote the respective components of the vectors. By squaring both sides we can also get rid of the square root: s^2 \cdot t^2 = (o_x + t \cdot v_x)^2+(o_y + t \cdot v_y)^2. Now we finally have a chance to sort out all the variables by simple transformations. The brackets can be squared out either manually or by binomic equivalences and then we can factor out t and order the terms to finally get t^2\cdot(v_x^2 + v_y^2 - s^2) + t\cdot2\cdot(o_x \cdot v_x + o_y \cdot v_y) + o_x^2 + o_y^2= 0 which is a quadratic equation that can be solved. Before doing so however we take a look at its components and try to simplify our calculation. Since we could fill in numbers for all variables but t we do so by substituting with simpler “helper” variables denoted by the following definitions: h_1 := v_x^2 + v_y^2 - s^2;h_2 := o_xv_x + o_yv_y;h_3 = o_x^2+o_y^2. Thus, by dividing our last equation by h_1 we get the simple quadratic equation t^2 +t\cdot 2\cdot h_2 + h_3 which, being quadratic, gives us two (or one, or none) solutions: t_{1,2} = -\frac{h_2}{h_1} \pm \sqrt{(\frac{h_2}{h_1})^2 - \frac{h_3}{h_1}}.

If the expression under the square root is negative the equation got no solution and we can never hit the object with our projectile since it is too slow. If the square root works out to some positive value we got two solutions. Solutions with a negative value for t can be thrown away too, unless you shoot with Tachyon-bullets 😉

There remains one case to be taken care of: h_1 might be 0 which would result in a division by 0 in the aproach shown above. But actually the case where h_1 is 0 is the mathematically simpler case. The quadratic equation collapses and all that remains is the simple equation t\cdot2\cdot h_2 + h_3 = 0 which of course results in t = -\frac{h_3}{2h_2}. Thanks to David (see comments below) for pointing that out. It got rid of an ugly ‘quick fix’ I used before.

And now in the form of Java code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
	/**
	 * Calculates the point of interception for one object starting at point
	 * <code>a</code> with speed vector <code>v</code> and another object
	 * starting at point <code>b</code> with a speed of <code>s</code>.
	 * 
	 * @see <a
	 *      href="http://jaran.de/goodbits/2011/07/17/calculating-an-intercept-course-to-a-target-with-constant-direction-and-velocity-in-a-2-dimensional-plane/">Calculating
	 *      an intercept course to a target with constant direction and velocity
	 *      (in a 2-dimensional plane)</a>
	 * 
	 * @param a
	 *            start vector of the object to be intercepted
	 * @param v
	 *            speed vector of the object to be intercepted
	 * @param b
	 *            start vector of the intercepting object
	 * @param s
	 *            speed of the intercepting object
	 * @return vector of interception or <code>null</code> if object cannot be
	 *         intercepted or calculation fails
	 * 
	 * @author Jens Seiler
	 */
	public static Point2D calculateInterceptionPoint(final Point2D a, final Point2D v, final Point2D b, final double s) {
		final double ox = a.getX() - b.getX();
		final double oy = a.getY() - b.getY();
 
		final double h1 = v.getX() * v.getX() + v.getY() * v.getY() - s * s;
		final double h2 = ox * v.getX() + oy * v.getY();
		double t;
		if (h1 == 0) { // problem collapses into a simple linear equation 
			t = -(ox * ox + oy * oy) / (2*h2);
		} else { // solve the quadratic equation
			final double minusPHalf = -h2 / h1;
 
			final double discriminant = minusPHalf * minusPHalf - (ox * ox + oy * oy) / h1; // term in brackets is h3
			if (discriminant < 0) { // no (real) solution then...
				return null;
			}
 
			final double root = Math.sqrt(discriminant);
 
			final double t1 = minusPHalf + root;
			final double t2 = minusPHalf - root;
 
			final double tMin = Math.min(t1, t2);
			final double tMax = Math.max(t1, t2);
 
			t = tMin > 0 ? tMin : tMax; // get the smaller of the two times, unless it's negative
			if (t < 0) { // we don't want a solution in the past
				return null;
			}
		}
 
		// calculate the point of interception using the found intercept time and return it
		return new Point2D.Double(a.getX() + t * v.getX(), a.getY() + t * v.getY());
	}
Leave a comment ?

21 Comments.

  1. Thanks, Jaran!

    This is exactly what I was looking for – you explain the maths really well (and I can see how it can be extended to 3D relatively simply).

    It’s helped me just code it up in C++.

    David

    • Thanks for your feedback. If you happen to write some piece of code that extends the maths to the 3rd dimension I’d appreciate if you’d posted it here 🙂

  2. I started expanding the code to work in 3D dimensions, but my problem was essentially still in a 2D plane so I ended up just using the two relevant components instead. Bits like the magnitude of the vector work out simply but I’m not sure about the calculation of the substitutions yet, hopefully there’s still a way to do it!

    One thing that caught me out is the case where ‘h1’ is zero. You adjust it to be a small number yourself in the code above, but the collision equation actually simplifies because the t^2 term disappears. You just get the linear equation:

    t*2*(ox*vx + oy*vy) + ox^2 + oy^2 = 0
    or
    t = -(ox^2 + oy^2) / (2 * (ox*vx + oy*vy))

    So that may save some processing in some cases – or seems a bit neater at any rate. 🙂

    • Hey David. Thanks for pointing that out. Of course you are right and I updated my article and the code. Much nicer now! 😀

      • Jaran,
        You have an order of operations bug in your above code for this fix. You need to add parenthesis around 2*h2. Line 32 should be:

        t = -(ox * ox + oy * oy) / (2*h2);

  3. niel jonny

    could some one explain what the small “a” and “b” mean and whats their purpose

  4. “a” and “b” are the position vectors that point from the origin of your coordinate system to the location of the object to be intercepted (A) and your “interceptor gun” (B).

  5. Joseph Jones

    That will work but you’re doing it the long way. There’s no need to find time and it really boils down to 1 simple equation (well 2). Here you are:

    y= mx+b : Definition of a line

    We need to find the same x when y’s are the same so set the equation up for x:
    m(A)x + b(A) = m(B)x + b(B)

    Solve for X:
    m(A)x – m(B)x = b(B) – b(A)
    x(m(A) – m(B)) = b(B) – b(A)
    x = (b(B) – b(A)) / (m(A) – m(B))

    We can find b:
    y=mx+b
    b=-mx+y

    Therefore by substitution:
    x = (-m(B)*(B)x + B(y) – (-m(A)*(A)x + A(y))) / (m(A) – m(B))

    We have all those variables so just solve for x.

    After we have x, just plug it into y=mx+b and solve for y.

    • In your line “x = …” what do you mean by (B)x and (A)x?
      And how do we have all those variables?

      You’re calculating the intersection of two geometric lines but are assuming that we already do have both m(a) and m(b)?

      But it is exactly the direction we have to point the gun at that we want to calculate and don’t have as a given.

  6. Joseph Jones

    Found the 2d, and here is some pseudocode for 3d (well it works in my game atleast, maybe it’ll help anyone that finds this forum):

    –P1 = Gun Position [x,y,z]
    –V1 = Gun Velocity [x,y,z]
    –P2 = Target Position [x,y,z]
    –V2 = Target Velocity [x,y,z]
    function Predict(P1,V1,P2,V2)
    AverageDistance = Distance(P1,P2)
    AverageSpeed = Speed(V2-V1)
    Time = AverageDistance /AverageSpeed
    Bullseye = P2 + V2*Time

    return Bullseye
    end

    function Speed(Velocity)
    return math.sqrt(Velocity.x^2 + Velocity.y^2 + Velocity.z^2)
    end

    function Distance(VectorA,VectorB)
    return math.sqrt((VectorA.x – VectorB.x)^2 + (VectorA.y – VectorB.y)^2 + (VectorA.z – VectorB.z)^2)
    end

    • When you subtract the speed vectors how do you know the vector of the bullet which is the unknown thing that we want to calculate? I can see how you can use both (scalar) speeds as an approximation though which, if done in every iteration of your game-loop might result in an acceptable result.
      So I guess as a first approximation you select a speed vector which points directly at the current position of the target?

  7. Hi

    I need some more help with the h1 = 0 scenario(s) if this discussion is still active:

    The case where target speed and bullet speed are both 0 is unhandled in your code. Not very interesting I guess, just pointing out.

    More interesting is the case where h1 is zero because target speed == bullet speed, and they are non-zero. This will occur *regardless* of the target vector direction. As far as I understand, t = -h3 / (2h2) assumes the target vector v to be pointing toward B. In this single case, the formula produces the correct result for t. For any other vector v, where |v| == s, this fails.

    • Hello Jan,

      I cannot see where t = -h_3 / (2h_2) makes an assumption about the direction of v towards B. This value for t is just, as you point out, the special result for h_1 := v_x^2 + v_y^2 - s^2 = 0 and can be derived because it reduces the first term (the quadratic term) in t^2\cdot(v_x^2 + v_y^2 - s^2) + t\cdot2\cdot(o_x \cdot v_x + o_y \cdot v_y) + o_x^2 + o_y^2= 0 to 0.

  8. You say “Thus, by dividing our last equation by h_1 we get the simple quadratic equation t^2 +t\cdot 2\cdot h_2 + h_3…”

    At this point, shouldn’t the equation be (using the notation above) t^2\cdot h_1 + t\cdot 2\cdot h_2 + h_3?

    This is further supported by the simplified equation (when h_1 = 0) of t\cdot 2\cdot h_2 + h_3 = 0, which removes the t^2\cdot h_1 because h_1 = 0.

    Apart from that, nice! 🙂 And very useful!

  9. Hi!
    I’m just here for the math, not the code, but can someone tell me if I’m missing something obvious here? Wouldn’t you need t to begin with to calculate Vx and Vy? This whole thing seems to assume that they are known, but if your target is moving, then they wouldn’t they also require t, making it unsolvable?

  10. Playing with Homing Missiles – Orbital BBQ Laser - pingback on 2017/07/21 at 19:59
  11. Exactly what I was looking for. Two points:

    – For the simple case, it can also occur that t < 0, so the test on t < 0 needs to be moved out of the else block.

    – Also, for the simple case, there is h2 in the denominator, which may be zero as well. An extra test for this would be in order.

    Otherwise, thanks, this is helping me out. The website url explains why.

  12. DoubleThought

    It is possible for both h1 and h2 to be zero, thus resulting in an uncaught division-by-zero error. However, if both are zero, then it is not possible for the bullet to catch up with the targeted object in any case (they will be the same speed, but the targeted object will be traveling perpendicular to the line between the starting point of the bullet and the starting point of the targeted object).

    I suggest inserting the below code before line 32:

    if (h2 == 0) {
    return null;
    }

  13. Just to confirm a few things,

    I agree with DoubleThought that h2 should be tested for 0.

    I agree with Mark that there is a typo in the explanation. It should say something like this: We get h1 times t squared plus 2 times t times h2 plus h3; and dividing that by h1 we get: t squared plus 2 times h2 over h1 plus h3 over h1.

    One new thing to contribute: I wanted to find the best intercept course that brings the bullet to within q units of distance of the moving target. This is what I came up with:

    var h1 = vx * vx + vy * vy – s * s
    var h2 = ox * vx + oy * vy – q * s
    var h3 = ox * ox + oy * oy – q * q

    After making that change to how h2 and h3 are calculated, it seems that the remaining code can remain the same.

  14. Really well explained. Thanks Jaran!
    Now how to take into account the gravitational force acting on y axis of B (the projectile)? Assume A (target) is unaffected by gravity.
    I am unable to work out the math, where to put the -1/2gt^2 part.

Reply to niel jonny ¬
Cancel reply


NOTE - You can use these HTML tags and attributes:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Trackbacks and Pingbacks: