Intercept an accelerating object with a constant velocity projectile

Question

I want to calculate the time that it takes to intercept a moving target "plane" with a projectile.

For example I have two 3 dimensional vectors (2 locations in 3D space).

Vector3 Person = {100, 0, 50};
Vector3 Plane = {200, 50, 100};

Now let's suppose the person has a gun that fires a projectile at 140M/s in a straight line with no acceleration or bullet drop.

I want to hit the plane, which has an initial velocity and a constant acceleration.

Vector3 VelocityPlane = {-10, -10, 0};
Vector3 AccelerationPlane = {0, -10, 0};

So what do I want to do:

Using the above information, I want to calculate where I should aim to hit the Plane.

To achieve that I found the following kinematic equation:

$$\text{Future Plane Position} = \text{Position} + \text{Velocity} \cdot \text {Time} + \left(\frac 1 2 \cdot \text{Acceleration} \cdot \text{Time}^2 \right)$$

In C++ Code it would look like this:

Vector calculatePositionToAimAt(double time, Vector startPosition, Vector velocity, Vector accel = 0) {
    return startPosition + (velocity * time) + (0.5 * accel * pow(time, 2));
}

The Problem:

This function above is great! But there is something subtle missing. I don't have the time to calculate the Position where I have to aim for the function above, which happens to be the most important piece of this problem.

For many, the easiest way would probably be simple: Time = Distance / Bulletspeed But that wont work in my case since the Plane has Acceleration and is moving.

So for that I found the following equation:

$$ \vec p_{\text{Plane}} + \vec v_{\text{Plane}} \cdot t + \left(\frac 1 2 \vec a \cdot t^2\right) = \vec p_{\text{Player}} + \vec v_{\text{Bullet}} \cdot t$$

Together with the explanation: Because time is not a vector and is a scalar, You no longer have to use vectors in this equation. You can just input the magnitudes of the vectors instead - which would look like this:

0.9 Seconds is the result. So this works fine.

Now the question is, how can I rewrite this in code so that it automatically solves such equations?

I already have found code that should solve this which looks like so:

double asymetricSingleAcceleratedInterceptionTime(Vector position1, Vector position2, Vector velocity1, Vector velocity2, Vector acceleration2) {
    /*
    sqrt( (2 * a2 * x1) + (2 * a2 * x2) + v1^2 - (2 * v1 * v2) + v2^2 ) - v1 + v2
    t = ------------------------------------------------------------------------------
    a2
    */

    if (acceleration2.x == 0 && acceleration2.y == 0 && acceleration2.z == 0) {
        return asymetricInterceptionTime(position1, position2, velocity1, velocity2);
    }

    return (sqrt(
        (2 * acceleration2 * position1).Length()
        + (2 * acceleration2 * position2).Length()
        + pow(velocity1.Length(), 2)
        - (2 * velocity1 * velocity2).Length()
        + pow(velocity2.Length(), 2))
        - velocity1.Length()
        + velocity2.Length()
        ) / (acceleration2).Length();
}

However, this always gave me wrong results, so I rewrote it a bit and this is what it looks like now:

Vector3 Player = Vector3(100, 0, 50);
Vector3 Plane = Vector3(200, 50, 100);
Vector3 BulletVelocit = Vector3(140,0,0);
Vector3 PlaneVelocit = Vector3(-10,-10,0);
Vector3 PlaneAcceleration = Vector3(0, -10, 0);

float pos1 = Player.Length();
float pos2 = Plane.Length();
float vel1 = BulletVelocit.Length();
float vel2 = PlaneVelocit.Length();
float acc1 = PlaneAcceleration.Length();


float calc1 = (2 * acc1 * pos1);
float calc2 = (2 * acc1 * pos2);
float calc3 = pow(vel1, 2);
float calc4 = (2 * vel1 * vel2);
float calc5 = pow(vel2, 2);
float calc6 = vel1;
float calc7 = vel2;
float calc8 = (acc1);

float time = (sqrt(calc1 - calc2 + calc3 - calc4 + calc5) + calc6 - calc7) / calc8;

std::cout << time;

This actually works but still has 2 problem:

1. it gives me wrong results when the velocity is minus, and I already know why
2. it gives me of course only one result and not 2, like above in the picture {0.9, 26.072}

About the first Problem: If I change this code from this:

float calc1 = (2 * acc1 * pos1); //acc1(acceleration) is +10 because of magnitude
float calc2 = (2 * acc1 * pos2);

to this:

float calc1 = (2 * -10 * pos1); //manually setting acceleration to -10 like in the example
float calc2 = (2 * -10 * pos2);

It gives me one right result: -26.072

So the second Problem.. I dont need the -26.072, but I need the other one which is 0.9

So my final questions are:

1. Is there other code that works better to give me the interception time based on velocity, acceleration etc like described above?

2. If not, how do I have to rewrite my code that it gives me the right positive result and not the negative one

3. How do I fix the magnitude acceleration problem which gives me wrong results.

Answer 1

This is exactly equivalent to this previous Q&A: Projectile Aim Prediction with Acceleration.

The only difference is that since your projectile doesn't accelerate, your \$\vec a_\text{Bullet} \$ term is \$\vec 0\$.

That gives us the following formula for the time-to-intercept \$t\$:

$$ \begin{align} 0 &= t^4 \left(\frac 1 4 \vec a \cdot \vec a \right)\\ &+ t^3 \left(\vec a \cdot \vec v \right)\\ &+ t^2 \left(\vec a \cdot \vec p + \vec v \cdot \vec v - s^2 \right)\\ &+ t^1 \left(2\vec p \cdot \vec v \right)\\ &+ 1 \left( \vec p \cdot \vec p\right)\\ \end{align}$$

Where \$\vec v, \vec a\$ are the velocity and acceleration of the plane, respectively, \$\vec p = \vec p_\text{Plane} - \vec p_\text{Player}\$ is the initial displacement of the plane from the player's firing position, and \$s\$ is the speed of your projectile.

Plug this into a quartic equation solver to get any positive time value solutions \$t_*\$, then plug each one in here to get the displacement from the firing position at the moment of interception:

$$\vec p_\text{interecept} = \vec p + \vec v t_* + \frac 1 2 \vec a t_*^2$$

And then use that to get your bullet velocity:

$$\vec v_\text{Bullet} = \frac 1 {t_*} \vec p_\text{interecept}$$