To make things a little simpler, let's work within the inertial frame of our launcher. That means our projectile will be firing from the point (0, 0, 0).
(You can subtract any arbitrary firing position from the target position to get this state of affairs for the purpose of the calculations, then add it back when you need to work with the resulting path)
v_p the launch velocity of our projectile relative to the launcher (vector)
t the time after launch when our projectile impacts the target (scalar)
s the speed of our projectile (scalar)
p_T the offset from the launcher to the target at the time of firing (vector)
v_T the velocity of the target relative to the launcher at the time of firing (vector)
a_T the constant acceleration of the target (vector)
a_p the constant acceleration of the projectile (vector)
Using our parametric equation for motion over time under constant acceleration...
p(t) = p(0) + v * t + (a/2.0) * t^2
We know that at the time of impact we want:
v_p * t + (a_p/2.0) * t^2 = p_T + v_T * t + (a_T/2.0) * t^2
Rearranging to isolate
v_p, we get:
v_p = (1/t)(t^2(a_T - a_p)/2.0 + t * v_T + p_T)
a = a_T - a_p since they'll appear paired from here on in.
Applying our speed constraint: (Here · refers to a dot product between vectors)
s = ||v_p||
s^2 = v_p·v_p
s^2 = (1/t^2)(t^4(a·a/4.0) + t^3(a·v_t) + t^2(a·p_T + v_T·v_T) + t(2.0 * v_T·p_T) + p_T·p_T)
t^4 * (a·a/4.0)
+ t^3 * (a·v_T)
+ t^2 * (a·p_T + v_T·v_T - s^2)
+ t * (2.0 * v_T·p_T)
That's a quartic equation in
t, which is not the most fun to solve. When the target is stationary, the
t terms vanish and we can use quadratic formula to solve for
t^2, but for the general case of a moving target there's not much we can do but bite the bullet.
You can either plug this into a math library to solve the polynomial, write your own polynomial solver, or approximate a solution using an iterative method like Newton's.
However you choose to go about it, you'll come out with up to 4 time values at which you can successfully intercept the target. Some could be negative, corresponding to projectiles being lobbed from the target sometime in the past to arrive at the firing position at time zero travelling at the desired speed - discard those and consider only the positive time values. (You might get none, if the target is simply out of firing range or moving too fast to intercept)
To turn any of these remaining valid time values into a launch velocity, do the following:
t value into the formula for the target's position
p(t) = p_T + v_T * t + (a_T/2.0) * t^2
This gives you the position of the target at the time of impact. Since this is the net offset in space that your projectile will have to bridge, we'll call this vector
Calculate the initial velocity (aim direction & speed) of an arc travelling to this point in space & time:
v_p = travel / t - (a_p/2.0) * t
By construction, launching the projectile with this velocity
v_p will intersect the target's movement while respecting the initial speed constraint we started with, to within the numerical accuracy of the operations we took to get here. (Some rounding can be introduced at each step)
Here's an animation of this formula at work:
This little test rig generates random target positions/velocities/accelerations and projectile gravity values, then calculates up to four intercept trajectories.
The white ball is the target, moving along its trajectory before t=0 (black curve) to its position at the time of firing (grey dot) and continuing along its future trajectory t > 0 (grey curve). Its path is marked by red dots at the solution points.
The intercepting projectiles and their trajectories appear in cyan & green at t = 0 and arc to intercept the target at their corresponding impact times.
For completeness, I also show the negative solutions if there are any, in yellow or magenta, firing from the target when it crosses a negative solution time and arriving at the firing point in the center at t = 0.