Problem

This is similar to a Closest Point of Approach calculation, but instead of the closest point(s), I want the points at which the objects are a given distance away from each other. (Imagine two spaceships moving on a collision course towards each other; I want to know at what point Ship B is within Ship A's weapons range.)

Attempted solution

I started by using this page which details a method of calculating the CPA. This seemed a promising start:

d(t) = |P(t) – Q(t)| = |w(t)| where w(t)=w0+t(u-v) with w0=P0-Q0

Rationale: for a desired distance of 5 units, I want to do:

5 = |w(t)|

Square, substitute variables, and expand. I'm working in 2D, so this isn't too bad. Variable substitution as follows:

w.x = a
w.y = b
t = x   //(Wolfram won't solve for other variables?!)
u.x = p
u.y = q
v.x = r
v.y = s


(a +x(p-r))^2 + (b + x(q-s))^2 = 25

Throw this into Wolfram Alpha and it gives the enormous formulas to solve for x.

And turn it into Python code:

def CPA_distance(e1, e2, dist):
"""Calculates the (two) times when e1 and e2 are exactly dist apart from each other."""
posd = e1.pos - e2.pos
dist2 = dist*dist

a = posd.x
b = posd.y
p = e1.vel.x
q = e1.vel.y
r = e2.vel.x
s = e2.vel.y

#TODO: account for the case where (p**2 - 2*p*r + q**2 - 2*q*s + r**2 + s**2) == 0 -> they'll never intersect

#NOTE: This function will return incorrect results if the velocities for the two entities are the same. In this case, simply check the distance between them.

denominator = (p**2 - 2*p*r + q**2 - 2*q*s + r**2 + s**2)
sqrt_part = sqrt(   (2*a*p - 2*a*r + 2*b*q - 2*b*s)**2
- 4*(a**2 + b**2 - dist2)*(p**2 - 2*p*r + q**2 - 2*q*s + r**2 + s**2))
sub = -a*p + a*r - b*q + b*s
t1 = ( -0.5 * sqrt_part + sub ) / denominator
t2 = (  0.5 * sqrt_part + sub ) / denominator

return (t1, t2)


This works in the majority of cases, but not all of them. It's possible for the calculation of sqrt_part to fail attempting to take the square root of a negative number, for instance with the following parameters (pseudocode):

e1 = ( pos = (5,5), vel = (0.5, 0.5) )
e2 = ( pos = (0,0), vel = (-0.5, 0.5) )
dist = 5


Question

What's going on that causes me to have a negative value inside the square root? Just drawing out the initial locations of e1 and e2 shows me that they'll never come within the required distance of each other but I'd have thought that I'd either get a negative time or get a divide by zero somewhere.

Is this a "valid" outcome which I can take to mean "will never get within distance, even in the past" or are my formulas all wrong?

• Just a note: the code would be simpler if you did veld = e1.vel - e2.vel and used that variable everywhere instead of p, r, q and s. – sam hocevar Mar 22 '14 at 11:53
• Agreed. However, Wolfram interprets ab as a*b, so I had to stick with single letter variables. – Floomi Mar 22 '14 at 13:45