Find the direction in which an object with linear motion should be thrown in order to touch with another object moving linearly

Question

Suppose we have two objects in a 3D space:

1. The first object, A, is at the given position P1 moving at a constant 3D Vector speed of V1
2. The second object, B, is at a given position that we will call P2

Knowing that B can only move at a speed of magnitude M1, in which direction should B move in order to touch A? The result must be expressed as a Unit 3D Vector

Here's a visual image I put up in Paint to help visualize the problem: Visual Example

Thank you.

Answer 1

Let the position of the target be a vector-valued function $$\vec{p}(t) = \vec{b} + t \vec{v} \tag{1a}\label{1a}$$ and the position of the interceptor $$\vec{P}(t) = \vec{B} + t \vec{V}, \quad \lVert \vec{V} \rVert = M \tag{1b}\label{1b}$$ so that at some time $\tau \gt 0$, the two collide: $$\vec{p}(\tau) = \vec{P}(\tau) \tag{2}\label{2}$$ Expanding and simplifying $\eqref{2}$ we get $$\begin{aligned} \vec{b} + \tau \vec{v} &= \vec{B} + \tau \vec{V} \\ \tau \vec{V} &= \vec{b} - \vec{B} + \tau \vec{v} \\ \vec{V} &= \vec{v} + \frac{\vec{b} - \vec{B}}{\tau} \\ \end{aligned} \tag{3}\label{3}$$ In Cartesian component form, we have four equations and four unknowns ($V_x$, $V_y$, $V_z$, and $\tau$): $$\left\lbrace ~ \begin{aligned} M^2 &= V_x^2 + V_y^2 + V_z^2 \\ V_x &= v_x + \frac{b_x - B_x}{\tau} \\ V_y &= v_y + \frac{b_y - B_y}{\tau} \\ V_z &= v_z + \frac{b_z - B_z}{\tau} \\ \end{aligned} \right. \tag{4a}\label{4a}$$ This can have zero, one, or two sets of real solutions: $$\begin{cases} \tau = \frac{\vec{v} \cdot (\vec{B} - \vec{b})}{\lVert\vec{v}\rVert^2 - M^2} \pm \frac{\sqrt{ M^2 \lVert \vec{B} - \vec{b} \rVert^2 - \lVert \vec{v}\times(\vec{B} - \vec{b})\rVert^2}}{\lVert \vec{v} \rVert^2 - M^2}, & \lVert\vec{v}\rVert \ne M \\ \tau = \frac{\lVert \vec{B} - \vec{b} \rVert^2}{2 \vec{v} \cdot (\vec{B} - \vec{b}) }, & \lVert\vec{v}\rVert = M \\ \end{cases} \tag{4b}\label{4b}$$ where only $\tau \gt 0$ makes sense. Note that the second only has a positive real solution when $\vec{v}$ is (in the half-space) towards $\vec{B}$ from $\vec{b}$. Then, $$\vec{V} = \frac{\vec{v} + \vec{b} - \vec{B}}{M \tau} \tag{4c}\label{4c}$$ and $\lVert\vec{V}\rVert = 1$ within rounding error.

You can obtain these using Maxima or wxMaxima (free, available for all operating systems) or some other Computer Algebra System like Maple, via e.g.

solve([ V_x^2 + V_y^2 + V_z^2 = M^2,
        V_x = v_x + (b_x - B_x)/tau,
        V_y = v_y + (b_y - B_y)/tau,
        V_z = v_z + (b_z - B_z)/tau ], [ tau, V_x, V_y, V_z ]);
solve([ V_x^2 + V_y^2 + V_z^2 = v_x^2 + v_y^2 + v_z^2,
        V_x = v_x + (b_x - B_x)/tau,
        V_y = v_y + (b_y - B_y)/tau,
        V_z = v_z + (b_z - B_z)/tau ], [ tau, V_x, V_y, V_z ]);

Here is an example Python 3 implementation of the above:

from math import sqrt

def intercept(m, dx, dy, dz, vx, vy, vz):
    tdiv = m*m - vx*vx - vy*vy - vz*vz
    while True:
        if tdiv != 0:
            tsqr = m*m*(dx*dx + dy*dy + dz*dz) - (dy*vz-dz*vy)**2 - (dz*vx-dx*vz)**2 - (dx*vy-dy*vx)**2
            if tsqr < 0:
                break
            tdif = sqrt(tsqr) / tdiv
            tsum = (vx*dx + vy*dy + vz*dz) / tdiv
            t1 = tsum + tdif
            t2 = tsum - tdif
            if t1 > 0 and t2 > 0:
                t = min(t1, t2)
            elif t1 > 0:
                t = t1
            elif t2 > 0:
                t = t2
            else:
                break
        else:
            tdiv = vx*dx + vy*dy + vz*dz
            if tdiv == 0:
                break
            t = 0.5 * (dx*dx + dy*dy + dz*dz) / tdiv
            if t <= 0:
                break

        return (dx+t*vx)/(t*m), (dy+t*vy)/(t*m), (dz+t*vz)/(t*m), t

    return None, None, None, None


def V(x, y, z):
    return "(%.6f, %.6f, %.6f)" % (x, y, z)

if __name__ == '__main__':
    from random import Random
    rng = Random()

    # Interceptor launch site
    origin_x = rng.uniform(-10, +10)
    origin_y = rng.uniform(-10, +10)
    origin_z = rng.uniform(-10, +10)

    # Interceptor speed
    speed = rng.uniform(2.0, 5.0)

    # Target starting location
    target_x0 = rng.uniform(-10, +10)
    target_y0 = rng.uniform(-10, +10)
    target_z0 = rng.uniform(-10, +10)

    # Target direction and speed
    target_vx = rng.uniform(-1, +1)
    target_vy = rng.uniform(-1, +1)
    target_vz = rng.uniform(-1, +1)
    target_v = sqrt(target_vx*target_vx + target_vy*target_vy + target_vz*target_vz)

    print("Interceptor launches from %s" % V(origin_x, origin_y, origin_z))
    print("               with speed %.6f" % speed)
    print("     Target launches from %s" % V(target_x0, target_y0, target_z0))
    print("     with velocity vector %s" % V(target_vx, target_vy, target_vz))
    print("                 at speed %.6f" % target_v)
    print("             in direction %s" % V(target_vx/target_v, target_vy/target_v, target_vz/target_v))

    dir_x, dir_y, dir_z, when = intercept(speed,
                                          target_x0 - origin_x,
                                          target_y0 - origin_y,
                                          target_z0 - origin_z,
                                          target_vx,
                                          target_vy,
                                          target_vz)
    if when is None:
        print("Interception is not possible.")
    else:
        print("   Interception occurs at %s" % V(target_x0 + when*target_vx,
                                                 target_y0 + when*target_vy,
                                                 target_z0 + when*target_vz))
        print("                        = %s" % V(origin_x + when*speed*dir_x,
                                                 origin_y + when*speed*dir_y,
                                                 origin_z + when*speed*dir_z))
        print("  if launced in direction %s" % V(dir_x, dir_y, dir_z))
        print("  whose Euclidean norm is %.9f" % sqrt(dir_x*dir_x + dir_y*dir_y + dir_z*dir_z))

When run, the above generates a random origin, random target, random target velocity, random interceptor speed, and finds the interceptor direction.