Projectile Aiming

The following code can be used to accurately determine the direction that a projectile should be fired at in order to collide with the target.

This code is capable of accurately targetting moving targets, falling targets, and projectiles uniformly accelerating in any direction. It was originally written for use in Deathball by UsAaR33

The Code

struct QuartSolution{        //solutions to quart!
  var float u[4];
};

//quadratic solver (using everybody's favorite algebra formula)
static function calculateQuadRoots(float a, float b, float c, out QuartSolution Q){
    local float sqterm;
    sqterm = b*b - 4*a*c;
    if (sqterm<0){ //imaginary root. return t=1
        Q.u[0]=-1;
        Q.u[1]=-1;
    }
    else{
        sqterm=sqrt(sqterm);
        a*=2;
        b*=-1;
        Q.u[0]=(b+sqterm)/(a);
        Q.u[1]=(b-sqterm)/(a);
    }
}

//Solve a 4th Degree Polynomial (Quartic) equation for 0s.
//taken from a javascript webpage (Explicitly states in source that source may be reused in any way)
//uses the quartic formula! :)
static function calculateQ(float aq, float bq, float cq, float dq, float eqin, out QuartSolution Q){
local float eq;
local float fq;
local float gq;
//local float hq;


// These are the squares of the local floatiables used to calculate the 4 roots-->
local float kq;
local float lq;
local float mq;
local float nq;
local float mq2;

local float compsw;
local float kqsw;
local float lqsw;
local float mqsw;

// Switch used in calculating REAL quartic roots)-->
local float sw;

// local floatiables for calculating REAL quartic roots)-->

local float kans;
local float lans;
local float theta;


local float x1;
local float x2;
local float x3;
local float x4;

local float x2a, x2b, x2c, x2d;
//local float x1b, x1b2, x2b2, x3b, x3b2, x4b, x4b2;
//more:
local float dnm;
local float a, b, c, d, f, g, h, k, m, m2,n,n2,r,rc;
local float calcy, calcp, calcr, calcq, calcx, calcmod;
local float dnmsw;
local int i;

// the 'q' suffix  denotes local floatiables used in the quartic equation
for (i=0;i<4;i++) 
    Q.u[i]=-1.0; //set to complex solutions
compsw=0;
kqsw=0;
lqsw=0;
mqsw=0;
dnmsw=0;
sw=0;


dnm=aq;      //note: this assumes aq is non-zero.  Of course it should be (eval 0.25g!)

//Simplifying by dividing all terms by the aq term called 'dnm' meaning denominator
aq=bq/dnm;
bq=cq/dnm;
cq=dq/dnm;
dq=eqin/dnm;
//Which yields an equation of the form X^4 + AX^3 + BX^2 + CX + D = 0

eq= bq-((3*aq*aq)/8);
fq= cq+ (aq*aq*aq/8) -(aq*bq/2);
gq= dq- (3*aq*aq*aq*aq/256) + (aq*aq*bq/16) - (aq*cq/4);

// SOLVING THE RESULTANT CUBIC EQUATION
// EVALUATING THE 'f'TERM

a=1; b=eq/2; c=((eq*eq)-(4*gq))/16; d= ((fq*fq)/64)*-1;

f = (((3*c)/a) - (((b*b)/(a*a))))/3;
//EVALUATING THE 'g'TERM

g = ((2*((b*b*b)/(a*a*a))-(9*b*c/(a*a)) + ((27*(d/a)))))/27;

//EVALUATING THE 'h'TERM
h = (((g*g)/4) + ((f*f*f)/27));
if (h > 0){
compsw=2;
m = (-(g/2)+ (sqrt(h)));
// K is used because math.pow cannot compute negative cube roots?
k=1;
if (m < 0) k=-1;
  else k=1;
m2 = ((m*k)**(1.0/3.0));

m2 = m2*k;
k=1;
n = (-(g/2)- (sqrt(h)));
if (n < 0)
  k=-1;
else
  k=1;

n2 = (n*k)**(1.0/3.0);
n2 *=k;
k=1;
kq=  ((m2 + n2) - (b/(3*a)));
kq=sqrt(kq);
// ((S+U)     - (b/(3*a)))
calcmod= sqrt((-1*(m2 + n2)/2 - (b/(3*a)))*(-1*(m2 + n2)/2 - (b/(3*a))) + (((m2 - n2)/2)*sqrt(3))*(((m2 - n2)/2)*sqrt(3)));
calcy=sqrt((calcmod-(-1*(m2 + n2)/2 - (b/(3*a))))/2);
calcx=(((m2 - n2)/2)*sqrt(3))/(2*calcy);
calcp=calcx+calcy;
calcq=calcx-calcy;
calcr=kq;

nq=(aq/4);
x1=kq+calcp+calcq-nq;
x4=kq-calcp-calcq-nq;


Q.u[0]=-x1; //appearently was incorrect by a factor of -1
Q.u[1]=-1; //complex
Q.u[2]=-1; //complex
Q.u[3]=-x4;
}


// FOR H < 0

if (h<=0){
r = sqrt((g*g/4)-h);
k=1;

if (r<0)
  k=-1;
// rc is the cube root of 'r'

rc = ((r*k)**(1.0/3.0))*k;
k=1;
theta =acos((-g/(2*r)));

kq= (2*(rc*cos(theta/3))-(b/(3*a)));

x2a=rc*-1;
x2b= cos(theta/3.0);
x2c= sqrt(3)*(sin(theta/3));
x2d= (b/3.0*a)*-1;

lq=(x2a*(x2b + x2c))-(b/(3*a));

mq=(x2a*(x2b - x2c))-(b/(3*a));

nq=(aq/4.0);
}

if (h<=0){

// psudo-fix 0 bug.. not the best.. but works
if (abs(kq)<1.0/(10000.0))
  kq=0;
if (abs(lq)<1.0/(10000.0))
  lq=0;
if (abs(mq)<-1.0/(10000.0))
  mq=0;
if (kq<0) {return;} else {kq=sqrt(kq);}
if (lq<0) {return;} else {lq=sqrt(lq);}
if (mq<0) {return;} else {mq=sqrt(mq);}

if (kq*lq>0){mq2=((fq*-1)/(8*kq*lq));kans=kq;lans=lq;}
if (kq*mq>0){mq2=((fq*-1)/(8*kq*mq));kans=kq;lans=mq;}
if (lq*mq>0){mq2=((fq*-1)/(8*lq*mq));kans=lq;lans=mq;}




if (compsw==0){
  x1=kans+lans+mq2-nq;
  Q.u[0]=x1;
  x2=kans-lans-mq2-nq;
  Q.u[1]=x2;
  x3=(kans*-1)+lans-mq2-nq;
  Q.u[2]=x3;
  x4=(kans*-1)-lans+mq2-nq;
  Q.u[3]=x4;
}
}
}

/*Calculate aiming ideal rotation for firing a projectile at a potentially moving target (assumes pawn physics)
 IN:
 -StartLoc = world location where projectile is starting at
 -EndLoc = world Location we wish to Target (should lie in the targetted actor)
 -ProjSpeed = speed of the projectile being fired
 -Gravity = a vector describing the gravity
 -Target = the actual targetted ACTOR
 -bLeadTarget = Can we track the target?  (the entire point of this function)
 OUT:
 -dest: Location where the projectile will collide with Target
 -returns vector describing direction for projectile to leave at
*/
static function vector GetShootVect(vector StartLoc, vector EndLoc, float ProjSpeed, vector Gravity, actor Target, bool bLeadTarget, out vector Dest)
{
  local QuartSolution Q;
  local float best, speed2D, HitTime;
  local vector Pr;
  local int i;
  local vector HitNorm, HitLoc;
    local vector D; //EndLoc-StartLoc
    local vector V; //Target.velocity

    D = EndLoc-StartLoc;
    V = Target.Velocity;
  //track falling actors
  if (bLeadTarget && Target.Physics==Phys_Falling){
            calculateQuadRoots(V dot V - ProjSpeed*ProjSpeed,
             2*(V dot D),D dot D,Q); //use quadratic formula
            for (i=0;i<2;i++)
            if (best<=0||(q.u[i]>0 && q.u[i]<best))
                best=q.u[i];
            Pr = normal(D/best + V)*ProjSpeed;
            if (best<=0 || Target.Trace(HitLoc,HitNorm,EndLoc+V*best+0.5*Gravity*best
             *best,EndLoc+vect(1,1,0)*V*best) == none){
              //will be falling:
                Dest = StartLoc + PR*best+0.5*Gravity*best*best;
              return normal(PR)*ProjSpeed;
            }
            else if (best>0)  //determine how long actor will be in air
                HitTime = vsize(HitLoc - (EndLoc+vect(1,1,0)*V*best))/vsize(vect(0,0,1)*V*best+0.5*Gravity*best);
            else
                HitTime = 0; //assume most time not in air?
  }

  //ASSUME GROUND TRACKING
    if (bLeadTarget && Target.Physics==Phys_Falling){   //trace down from target to get ground normal
    Target.Trace(HitLoc,HitNorm,EndLoc+normal(Gravity)*5000,EndLoc);
        D.z=HitLoc.z-StartLoc.Z;  //set destination.z to floor, wipe out velocity.z and re-eval assuming ground
        V.z=0;    //no longer falling - view velcocity in 2D
        if (HitTime>0.5){  //True if likely in air most of time (in which case keep current V.X and V.y)
            V.z -= HitNorm.Z * (V dot HitNorm);
        }
        else{ //otherwise alter all of velocity vector, but keep current 2D speed
            speed2D = vsize(V);
            V=normal(V)*speed2D; //assume the same x and y speed if in air most time
            V -= HitNorm * (V dot HitNorm);   //recalculate players velocity on a slope using hitnormal  (assumes v.x and v.y is "ground speed")
            V=normal(V)*speed2D; //assume the same x and y speed if in air most time
        }
  }
    //todo: add traces to check side walls?
    //note: walking velocity *should* factor in current slope
    best=0;
    if (bLeadTarget && V!=vect(0,0,0)){
        calculateQ(0.25*(Gravity dot Gravity),(-Gravity) dot V,(-Gravity) dot D +
          V dot V - ProjSpeed*ProjSpeed,2*(V dot D),D dot D,Q);
        for (i=0;i<4;i++)
        if (best<=0||(q.u[i]>0 && q.u[i]<best))
            best=q.u[i];
    }
    else{ //don't lead. assume stationary target
            calculateQuadRoots(0.25*(Gravity dot Gravity),(-Gravity) dot D - ProjSpeed*ProjSpeed,D dot D,Q);
            for (i=0;i<2;i++)
            if (best<=0||(q.u[i]>0 && q.u[i]<best))
                best=q.u[i];
      if (best>0)
                best=sqrt(best);
  }

  if (best<=0){   //projectile is out of range
      //Warning: Out of range adjustments assume gravity is parallel to the z axis and pointed downward!!
        Pr.z =ProjSpeed/sqrt(2.0); //determine z direction of firing
        best = -2*Pr.z/Gravity.z;
        best+=(vsize(D)-pr.Z*best)/ProjSpeed; //note p.z = 2D vsize(p)  (this assumes ball travels in a straight line after bounce)
        //now recalculate PR to handle velocity prediction (so ball at least partially moves in direction of player)
    Pr = D/best + V - 0.5*Gravity*best;
        //now force maximum height again:
        Pr.z=0;
        Pr = (ProjSpeed/sqrt(2.0))*normal(Pr);
        Pr.z = ProjSpeed/sqrt(2.0); //maxmimum
        Dest = StartLoc + PR*best+0.5*Gravity*best*best;
        return Pr;
  }

  Pr = normal(D/best + V - 0.5*Gravity*best)*ProjSpeed;

    Dest = StartLoc + PR*best+0.5*Gravity*best*best;
    return Pr;
}

How To Use It

Simply call the function GetShootVec() with the following parameters:

StartLoc is the world location where the projectile is going to be spawned at.

EndLoc is the world Location that is being targetted (it should lie within the targetted actor)

ProjSpeed is speed at which the projectile being fired at

Gravity is a vector describing any uniform acceleration the projectile is under (gravity, "english", etc.)

Target is the actor being targetted. Generally it should be a pawn, but any class will work.

bLeadTarget. If this is true, then "velocity prediction" will occur (that is the projectile will be fired ahead of the target and still collide withit, assuming its velocity does not change in a nonpredictable way. If false, the projectile is guarenteed to collide with the targetted location

The function will return a vector in the direction the projectile be launched at and with a magnitude equal to the speed of the projectile (as passed in the parameters).

In addition, a vector dest is returned. This specifies the world location where the projectile will collide with the Target (if the target is moving and bLeadTarget is enabled this vector will differ from endloc).

Notes

Speed

If the target is stationary or bLeadTarget is false, then the quadratic formula can be applied to the 4th degree polynomial that is the core of this algorithm. (the coefficients of t^3 and t^1 become 0). The algorithm takes roughly 0.20 ms to complete (on an AMD Athlon XP 2100+ test machine), equivilent to a drop from 80.0 to 78.7 fps if this GetShootVec() is called once per frame.

If velocity prediction is needed, then the quartic formula (calculateQ) must be applied. This increases the algorithm's time to 0.28 ms (on an AMD Athlon XP 2100+ test machine),, equivilent to a drop from 80.0 to 78.2 fps if this GetShootVec() is called once per frame.

Falling Targets

If a target is falling, the algorithm will verify that the target will not land before the projectile reaches it. If it is determined that the target will land, then the target's position will be re-evaluated, assuming any movement made on the ground is modeled under phys_walking. (thus, this algorithm likely will cause the fired projectile to miss a bBounce=true projectile target).

Please note that it is assumed that the net acceleration of the target is equal to the gravity vector of the projectile. If this is not the case, the code should not be too hard to modify.

Out of Range Targets

Should a target be out of range, this algorithm will attempt to bring the projectile as close as possible to the target. It does make the "appoximate" assumption that the projectile being fired will travel in an "approximately" straight line at its specifified speed toward the target.

Furthermore, the out-of-range routine assumes that the projectile's acceleration direction passed is parallel to the Z axis and pointed downward. It should not be too hard to extend this code to handle different acceleration.

Projectiles under constant acceleration

If you wish your projectile to undergo other constant acceleration (that is curving), this is entire possible. Simply be sure that the gravity vector reflects this. Please note though that the Out Of Range routine does not properly handle this.

That structure that has nothing more than an array

I was thinking of eventually adding other variables to that structure, but I'm not sure if that will ever be needed. Feel free to remove the structure altogether and just make the variable the array.

The Quartic Equation solver

It is converted javascript code that is available on the internet (used with permission). Note that the code is highly unoptimized. Feel free to fix that.

How does this work?

Physics

Basically nothing more than simple kinematics.

Abbriviation notes:

D=Displacement

Vp=launch velocity of Projectile

Vps = launch speed of projectile

Va=velocity of target

t=time

G=acceleration of projectile

The displacement of a projectile can be described with the following vector equation:

D=Vp*t + 0.5*G*t^2

The size of V must be the speed of the projectile.

From this, two equations (four in scalar form) can be derived:

vsize(Vp) = Vp

Vp*t + 0.5*G*t^2 = Va*t + D

(Note that in the case where the target is also falling (under the same acceleration as the projectile), G becomes 0.)

After a bunch of math, the following equation pops out:

0=0.25*(G dot G)*t^4+((-G) dot Va)*t^3 + ((Va dot Va) - Vps^2 + (-G) dot D) * t^2 + 2*(Va dot D)*t + (D^2)

This is a 4th degree polynomial, which can be solved with the Quartic equation.

In the case when Va = 0 (stationary or no leading), the t and t^3 terms become 0, allowing the much nicer quadratic equation to be used. (the result simply needs to be squarerooted).

Once the time is found, the launch velocity is easily calculated (by using the above simultaneous equations):

  Vp= normal(D/lowestpostitivetime + Vp - 0.5*Gravity*lowestpostitivetime )*Vps

And that's pretty much it. The details of the code should be explained well enough through comments.