I'm trying to determine direction in which a missile should accelerate to hit a target in a 3D space simulation. The position and velocity of the missile, as well as the position, velocity, and acceleration of the target are known. Additionally, the missile has a thrust that it can produce in any direction. I've come up with what, as far as I can tell, should be a working solution, but in code, it's not giving a valid result (more info later). Here are the variables I have defined:
$p_t$ - Current Target Position (3-vector, known)
$v_t$ - Current Target Velocity (3-vector, known)
$a_t$ - Current Target Acceleration (3-vector, known)
$p_m$ - Current Missile Position (3-vector, known)
$v_m$ - Current Missile Velocity (3-vector, known)
$a_m$ - Desired Missile Acceleration (3-vector, unknown)
$|a_m|$ - Magnitude of Missile Acceleration (known)
I set up these two equations based on the known information:
$$p_m + v_m * t + 0.5 * a_m * t^2 = p_t + v_t * t + 0.5 * a_t * t^2$$
$${{a_m}_x}^2+{{a_m}_y}^2+{{a_m}_z}^2=|a_m|^2$$
Now, I try to solve this system of equations to determine a time and acceleration direction where the missile could intersect with the target while accelerating at the given rate. First, solve the first equation for $a_m$:
$$0.5 * a_m * t^2 = p_t + v_t * t + 0.5 * a_t * t^2 - p_m - v_m * t$$
$$a_m = \frac{p_t + v_t * t + 0.5 * a_t * t^2 - p_m - v_m * t}{0.5 * t^2}$$
To simplify slightly, I define two constants for relative velocity and position:
$C_1=p_t-p_m$
$C_2=v_t-v_m$
$$a_m = \frac{C_1 + C_2 * t + 0.5 * a_t * t^2}{0.5 * t^2}$$
This equation can be split into x, y, z components to get three separate equations solved for ${{a_m}_x}$, ${{a_m}_y}$, and ${{a_m}_z}$. For each component $n$, I define a new variable $S_n$ as the numerator of this equation:
$$S_n = {C_1}_n + {C_2}_n * t + 0.5 * {a_t}_n * t^2$$
And consequently,
$${a_m}_n = \frac{S_n}{0.5 * t^2}$$
This makes it simpler to plug in to the other equation ${{a_m}_x}^2+{{a_m}_y}^2+{{a_m}_z}^2=|a_m|^2$:
$$(\frac{S_x}{0.5 * t^2})^2+(\frac{S_y}{0.5 * t^2})^2+(\frac{S_z}{0.5 * t^2})^2=|a_m|^2$$
$$\frac{{S_x}^2+{S_y}^2+{S_z}^2}{(0.5 * t^2)^2}=|a_m|^2$$
$${S_x}^2+{S_y}^2+{S_z}^2=0.25*|a_m|^2*t^4$$
I'm not sure if there's a way to proceed without expanding $S_n$ into a big mess, so that's what I'm going to do.
$$({C_1}_x + {C_2}_x * t + {a_t}_x * t^2 * 0.5)^2 +$$
$$({C_1}_y + {C_2}_y * t + {a_t}_y * t^2 * 0.5)^2 +$$
$$({C_1}_z + {C_2}_z * t + {a_t}_z * t^2 * 0.5)^2 = 0.25*|a_m|^2*t^4$$
And then multiplying out:
$$({C_1}_x^2 + {C_1}_x*{C_2}_x*t + {C_1}_x*{a_t}_x*t^2*0.5) +$$
$$({C_2}_x*t*{C_1}_x + {C_2}_x^2*t^2 + {C_2}_x*t^3*{a_t}_x*0.5) +$$
$$({a_t}_x*t^2*0.5*{C_1}_x + {a_t}_x*t^3*0.5*{C_2}_x + {a_t}_x^2*t^4*0.25) +$$
$$({C_1}_y^2 + {C_1}_y*{C_2}_y*t + {C_1}_y*{a_t}_y*t^2*0.5) +$$ $$({C_2}_y*t*{C_1}_y + {C_2}_y^2*t^2 + {C_2}_y*t^3*{a_t}_y*0.5) +$$ $$({a_t}_y*t^2*0.5*{C_1}_y + {a_t}_y*t^3*0.5*{C_2}_y + {a_t}_y^2*t^4*0.25) +$$
$$({C_1}_z^2 + {C_1}_z*{C_2}_z*t + {C_1}_z*{a_t}_z*t^2*0.5) +$$ $$({C_2}_z*t*{C_1}_z + {C_2}_z^2*t^2 + {C_2}_z*t^3*{a_t}_z*0.5) +$$ $$({a_t}_z*t^2*0.5*{C_1}_z + {a_t}_z*t^3*0.5*{C_2}_z + {a_t}_z^2*t^4*0.25) =$$
$$0.25*|a_m|^2*t^4$$ This can be reorganized to group corresponding terms by component: $$({C_1}_x^2 + {C_1}_y^2 + {C_1}_z^2) +$$ $$({C_1}_x*{C_2}_x*t + {C_1}_y*{C_2}_y*t + {C_1}_z*{C_2}_z*t) +$$ $$({C_1}_x*{a_t}_x*t^2*0.5 + {C_1}_y*{a_t}_y*t^2*0.5 + {C_1}_z*{a_t}_z*t^2*0.5) +$$
$$({C_2}_x*t*{C_1}_x + {C_2}_y*t*{C_1}_y + {C_2}_z*t*{C_1}_z) +$$ $$({C_2}_x^2*t^2 + {C_2}_y^2*t^2 + {C_2}_z^2*t^2) +$$ $$({C_2}_x*t^3*{a_t}_x*0.5 +{C_2}_y*t^3*{a_t}_y*0.5 + {C_2}_z*t^3*{a_t}_z*0.5) +$$
$$({a_t}_x*t^2*0.5*{C_1}_x + {a_t}_y*t^2*0.5*{C_1}_y + {a_t}_z*t^2*0.5*{C_1}_z) +$$ $$({a_t}_x*t^3*0.5*{C_2}_x + {a_t}_y*t^3*0.5*{C_2}_y + {a_t}_z*t^3*0.5*{C_2}_z) +$$ $$({a_t}_x^2*t^4*0.25 + {a_t}_y^2*t^4*0.25 + {a_t}_z^2*t^4*0.25) =$$
$$0.25*|a_m|^2*t^4$$ From here, each line of the equation can be re-written as a dot product: $$(C_1 \cdot C_1) + t(C_1 \cdot C_2) + 0.5*t^2(a_t \cdot C_1) +$$ $$t(C_2 \cdot C_1) + t^2(C_2 \cdot C_2) + t^3*0.5(C_2 \cdot a_t) +$$ $$t^2*0.5(a_t \cdot C_1) + t^3*0.5(a_t \cdot C_2) + t^4*0.25(a_t \cdot a_t) -$$
$$0.25*|a_m|^2*t^4 = 0$$ And the terms can be grouped based on the power of t, into a quartic equation of the form $a + bt + ct^2 + dt^3 + et^4$, where: $$a = C_1 \cdot C_1$$ $$b = 2(C_1 \cdot C_2)$$ $$c = (a_t \cdot C_1) + (C_2 \cdot C_2)$$ $$d = C_2 \cdot a_t$$ $$e = 0.25(a_t \cdot a_t) - 0.25*|a_m|^2$$ I can solve this quartic equation, and take the lowest positive, real root, if there is one, to be the time at which the missile intercepts the target. The time can then be used to calculate $a_m$ by plugging it back in. Ok, so, in theory, as far as I can tell, this should work correctly. However, when I actually try it in my simulation, the resulting acceleration is clearly invalid, as it doesn't come close to satisfying the equation ${{a_m}_x}^2+{{a_m}_y}^2+{{a_m}_z}^2=|a_m|^2$ (it will have a magnitude that's way smaller than the expected $|a_m|$. I've gone over my solution several times, and confirmed that my code matches. I do know that my quartic solver is not the problem, since it gets the same answer as Wolfram and Desmos when given the quartic equation shown above. Can anyone see what might be going wrong in my process here?
I actually figured this out. The problem was in the order of arguments to my quartic solver; the math seems to be correct. I'll leave this up in case anybody might find it useful in solving a similar problem.