# Analytical Ballistic Trajectories with Approximately Linear Drag

Authors: Giliam J. P. de Carpentier
Journal: International Journal of Computer Games Technology
 (2014)

---

## Abstract

This paper introduces a practical analytical approximation of projectile trajectories in 2D and 3D roughly based on a linear drag model and explores a variety of different planning algorithms for these trajectories. Although the trajectories are only approximate, they still capture many of the characteristics of a real projectile in free fall under the influence of an invariant wind, gravitational pull, and terminal velocity, while the required math for these trajectories and planners is still simple enough to efficiently run on almost all modern hardware devices. Together, these properties make the proposed approach particularly useful for real-time applications where accuracy and performance need to be carefully balanced, such as in computer games.


---

## Body

## 1. Introduction

A ballistic trajectory is the path of an object that is dropped, thrown, served, launched, or shot but has no active propulsion during its actual flight. Consequently, the trajectory is fully determined by a given initial velocity and the effects of gravity and air resistance. Mortars, bullets, particles, and jumping computer game characters (between key presses) are all examples of ballistics, while actively controlled aircraft and rocket-propelled grenades are not.Describing the exact motion of an object in free fall is a classic problem that can become quite complex when including effects like drag, turbulence, height-dependent medium pressure, position-dependent gravity, buoyancy, lift, and rotation. For this research paper, the problem will be approached on a relatively pragmatic level, as it will be based on a reasonably simple drag model that does not consider the dynamics of projectile rotation and assumes that wind, gravity, and terminal velocity all remain fixed over the whole trajectory. As a result, some accuracy will be sacrificed for the sake of computational efficiency and flexibility in practical use, while still maintaining much of the essence of the ballistic motion through a resistive medium. Although such a choice does not make much sense for most scientific and military applications, it does make sense for computer games, where performance is typically more important than physical correctness.Currently, computer games hardly ever use trajectories influenced by air resistance when complex planning is required. That might partially be because implementing a computer player that is capable of quickly calculating the ideal angle to fire a mortar with a fixed initial speed to hit a given target, for example, is harder when having to take into account drag and wind conditions. In fact, the added complexity and computational intensity that would be required for working with many of the current drag models might simply not be justifiable.This paper introduces a trajectory model that is designed to fit in the gap where working with accurate models would be too complex and working with simple dragless parabola-shaped trajectories would be insufficient. The proposed model’s use and practicality will be demonstrated by covering a number of its properties and showing how these trajectories can be planned in a variety of ways.In Section2, previous work is shortly considered. In Section 3, the new model will be introduced and qualitatively compared to other models. In Section 4, the first planner is covered. Other planners use a different space as explained in Section 5 and will be covered in Section 6. This is followed by a discussion of Future Work in Section 7, the Conclusion in Section 8, and the References. The included Algorithms 1 and 2 contain C++ functions that efficiently solve the two most complex planning problems.Algorithm 1:Specialized C++ Quartic Solver.
d
ouble GetTimeToTargetRWithMinimalInitialSpeed(double k, double vInfinity,
double rX, double rY){
// 1. Start by getting coefficients for the function f(t) = a4*t
∧
4 + a3*t
∧
3
// + a1*t + a0 which is 0 at the sought time-to-target t. Solving f(t) = 0
// for t > 0 is equivalent to solving e(u) = f(1/u)*u
∧
4 = a0*u
∧
4 + a1*u
∧
3 +
// a3*u + a4 = 0 for u where u = 1 / t, but the latter is more well-behaved,
// being a strictly concave function for u > 0 for any set of valid inputs,
// so solve e(u)=0 for u instead by converging from an upper bound towards
double kVInfinity = k * vInfinity, rr = rX * rX + rY * rY; /
/ the root and
double a0 =
−
rr, a1 = a0 * k, a3 = k * kVInfinity * rY; // return 1
/
u.
double a4 = kVInfinity * kVInfinity;
double maxInvRelError = 1.0E6; // Use an achievable inverse error bound.
double de, e, uDelta = 0;
// 2. Set u to an upper bound by solving e(u) with a3 = a1 = 0, clamped by
// the result of a Newton method’s iteration at u = 0 if positive.
double u = std::sqrt(kVInfinity / std::sqrt(rr));
if (rY < 0) u = std::min(u,
−
vInfinity / rY);
// 3. Let u monotonically converge to e(u)’s positive root using a modified
// Newton’s method that speeds up convergence for double roots, but is likely
// to overshoot eventually. Here,“e” =  e(u) and“de” = de(u)
/
du.
for (int it = 0; it < 10; ++it, uDelta = e / de, u
−
= 1.9
* uDelta){
de = a0 * u; e = de + a1; de = de + e; e = e * u;
de = de * u + e; e = e * u + a3; de = de * u + e; e = e * u + a4;
if (
!
(e < 0
&
& de < 0)) break;
// O
vershot the root.
}
u += 0.9 * uDelta; // Trace back to the unmodified Newton method’s output.
// 4. Continue to converge monotonically from the overestimated u to e(u)’s
// only positive root using Newton’s method.
for (int it = 0; uDelta * maxInvRelError > u && it < 10; ++it){
de = a0 * u; e = de + a1; de = de + e; e = e * u;
de = de * u + e; e = e * u + a3; de = de * u + e; e = e * u + a4;
uDelta = e / de; u
−
= uDelta;
}
// 5. Return the solved time t to hit [rX, rY], or 0 if no solution exists.
return u > 0
? 1 / u :  0;
}Algorithm 2:Specialized C++ Quartic Solver.
d
ouble GetTimeToTargetRGivenInitialSpeedS(double k, double vInfinity, double rX,
double rY, double s, bool highArc){
// 1. Start by getting coefficients for the function f(t) = a4*t
∧
4 + a3*t
∧
3
// + a2*t
∧
2 + a1*t + a0 which is 0 at the sought time-to-target t. Solving
// f(t) = 0 for t > 0 is equivalent to solving e(u) = f(1/u)*u
∧
3 = a0*u
∧
3 +
// a1*u
∧
2 + a2*u + a3 + a4/u for u where u = 1 / t, but the latter is more
// well-behaved, being a strictly convex function for u > 0 for any set of
// inputs iff a solution exists, so solve for e(u) = 0 instead by converging
// from a high or low bound towards the closest root and return 1/u.
double kRX = k * rX, kRY = k * rY, kRXSq = kRX * kRX, sS = s * s;
double twoKVInfinityRY = vInfinity * (kRY + kRY), kVInfinity = k * vInfinity;
double a0 = rX * rX + rY * rY, a1 = (k + k) * a0;
double a2 = kRXSq + kRY * kRY + twoKVInfinityRY − sS;
double a3 = twoKVInfinityRY * k, a4 = kVInfinity * kVInfinity;
double maxInvRelError = 1.0E6;   // Use an achievable inverse error bound.
double maxV0YSq = sS − kRXSq;// maxV0YSq is the max squared  “V0.y”  that leaves
double e, de, u, uDelta = 0;  // enough  “V0.x” to reach rX horizontally.
// 2. Set u to a lower/upper bound for the high/low arc, respectively.
if (highArc){// Get smallest u vertically moving rY at max possible +v0.y.
double minusB = std::sqrt(maxV0YSq) − kRY;
double determ = minusB * minusB − (twoKVInfinityRY + twoKVInfinityRY);
u = (kVInfinity + kVInfinity)
/ (minusB + std::sqrt(determ));
maxInvRelError = −maxInvRelError; // Convergence over negative slopes.
} else if (rY < 0){// Get largest u vertically moving rY at most neg. v0.y.
double minusB = −std::sqrt(maxV0YSq) − kRY;
double determ = minusB * minusB − (twoKVInfinityRY + twoKVInfinityRY);
u = (minusB − std::sqrt(determ)) / (rY + rY);
// Clamp the above bound by the largest u that reaches rX horizontally.
u = std::min(s / rX − k, u);
} else u = s / std::sqrt(a0) − k; // Get the (largest) u hitting rX
// horizontally a.s.a.p. while launching in the direction of [rX,rY].
// 3. Let u monotonically converge to e(u)'s closest root using a modified
// Newton's method, almost scaling the delta as if the solution is a double
int it = 0; // root. Note that  “e”  = e(u) * u
∧
2 and  “de”  = de(u)/du * u
∧
2.
for (; it < 12; ++it, uDelta = e / de, u −= 1.9 * uDelta){
de = a0 * u; e = de + a1; de = de + e; e = e * u + a2; de = de * u + e;
e = e * u + a3; e = (e * u + a4) * u; de = de * u * u − a4;
if (!(u > 0 && de * maxInvRelError > 0 && e > 0)) break; // Overshot.
}
u += 0.9 * uDelta; // Trace back to unmodified Newton method's output.
// 4. Continue to converge monotonically to e(u)'s closest root using
// Newton's method from the last known conservative estimate on the convex
// function. (Note that in practice, u will have converged enough in <12
for (; u > 0 && it < 12; ++it){
// iterations iff a solution does exists.)
de = a0 * u; e = de + a1; de = de + e; e = e * u + a2; de = de * u + e;
e = e * u + a3; e = (e * u + a4) * u; de = de * u * u − a4;
uDelta = e / de; u −= uDelta;
if (!(de * maxInvRelError > 0)) break; // Wrong side of the convex“dip”.
if (uDelta * maxInvRelError < u && u > 0) return 1 / u; // 5a. Found it!
}
// 5b. If no solution was found, return 0. This only happens if s (minus
// a small epsilon) is too small to have a solution, the target is at the
return 0; // origin, or the parameters are so extreme they cause overflows.
}


## 2. Previous Work

The motion of ballistic projectiles has been covered in many physics papers and textbooks, and all of these use their own set of assumptions to create an approximate model of the forces acting on a projectile. Although some ballistics research focuses on effects like shape and orientation [1], spin [2], or (sub)orbital flight [3], most works on ballistic trajectories assume a fixed gravitational pull and use a simplified drag model. These drag models typically ignore all effects of in-flight rotation and are only dependent on the local velocity relative to the medium and on a given fixed terminal velocity or drag coefficient.The used drag model influences both realism and computational complexity. For example, when no drag force is applied, the trajectory will always be a parabola, which is easy to work with and plan for. If the drag force is chosen to be linear in velocity, an explicit function describing the trajectory can be found by solving a set of linear differential equations [4]. This transcendental function is already computationally harder to calculate and even harder to plan with (i.e., solve for) [5]. To approximate reality even better, the drag force can be made quadratic in the object’s velocity relative to the medium. But as no exact analytic solution for the resulting trajectory exists, calculating a trajectory requires either crude approximation or numerical integration, and planning a trajectory requires reiteration [6–8].The research in this paper is based on a novel approximation of the trajectory function that follows from the linear drag model, sacrificing some of its moderate accuracy for a further increase in both efficiency and flexibility.


## 3. The Approximated Trajectory

### 3.1. Ballistic Parameters

Before presenting the proposed trajectory function and its properties, the necessary ballistic parameters will be defined here first. Starting with a general note, vector variables in this paper are always distinguished from scalar variables by the→ symbol above their names. Also, the length of any vector a
→ is denoted as ∥
a
→
∥, which is equal to a
→
·
a
→, where · is the dot product.A ballistic object is assumed to be launched from the initial positionp
→
0 with the initial velocity v
→
0 at time t
=
0. Furthermore, the object will be travelling through a medium (e.g., the air) which itself travels at the fixed (wind) velocity v
→
medium. It is also pulled by gravity at the fixed gravitational acceleration g, which is roughly 9.81 m/s2 for “earthly” applications. The amount of drag k while moving through the medium is defined as follows.
(1)
k
=
1
2
g
∥
v
→
terminal
∥
.Here,v
→
terminal is the invariant terminal velocity relative to v
→
medium that is reached eventually as the forces of gravity and air resistance finally cancel each other out. On earth, that is equivalent to saying that v
→
terminal is the fixed velocity that is approached when the object is dropped from an enormous height on a windless day. Lastly, the absolute terminal velocity v
→
∞, being the absolute velocity approached when time t goes to ∞, is therefore
(2)
v
→
∞
=
v
→
terminal
+
v
→
medium
.Together,p
→
0, v
→
0, v
→
∞, and k uniquely define a trajectory in the proposed model. To give a real-world example of the parameters defined above, suppose a tennis ball with a terminal velocity of 30 m/s is served at 50 m/s at a 30° angle from the 381 m high roof of the Empire State building into a 10 m/s horizontal wind in 2D. Then, p
→
0
=
[
0
,
381
], v
→
0
=
[
50
cos
⁡
(
3
0
°
)
,
50
sin
(
3
0
°
)
]
≈
[
43.30
,
25
], k
=
0.1635, and v
→
∞
=
[
-
10
,
-
30
]. The trajectory resulting from these values is shown in Figure 1.Figure 1
Serving a tennis ball from the Empire State Building upwind.


### 3.2. Deriving the Trajectory Function

In terms of the parameters defined above, the differential equation for the exact linear drag model has the following analytic solution:(3)
p
→
linear
(
t
)
=
1
2
k
(
v
→
0
-
v
→
∞
)
(
1
-
e
-
2
k
t
)
+
v
→
∞
t
+
p
→
0
.This function calculates the 2D or 3D positionp
→
linear on a trajectory at time t where t
≥
0. The above function is far from new and will not be explained here in detail, as it has already been covered in many textbooks [4], occasionally even targeting game developers in particular [9] (albeit with slightly different notation and parameter definitions).The function above will not be used directly in this paper. Instead, it will be approximated by substituting its exponential functione
x with the first degree rational function (
2
+
x
)
/
(
2
-
x
) shown in Figure 2. One of the reasons for selecting this approximation over all possible other approximations to e
x is that it has a value, first derivative and second derivative that match those of e
x at x
=
0. This means that it approximates e
x near x
=
0 well and therefore will guarantee a good approximation of p
→
linear
(
t
) near t
=
0. Furthermore, the first derivative of (
2
+
x
)
/
(
2
-
x
) monotonically decreases from 1 to 0 as x tends from 0 to -
∞, similar to the first derivative of e
x itself. When used to approximate e
x in p
→
linear
(
t
), this property will cause the initial velocity to monotonically converge to the terminal velocity over time. Note that no polynomial approximation of e
x has this specific property, and of all possible rational functions that do possess the above properties, the proposed approximation is the simplest and thus the most efficient. Lastly, its inverse is also a first degree rational function, resulting in relatively simple algebraic solutions for all (otherwise algebraic) equations that use it to approximate e
x.Figure 2
Approximatinge
x using a simple rational function.When the exponential function inp
→
linear
(
t
) is substituted by the rational approximation, the following function is the result
(4)
p
→
(
t
)
=
(
v
→
0
+
k
t
v
→
∞
)
t
1
+
k
t
+
p
→
0
.Because of the aforementioned properties,p
→
(
t
) will not only be more efficient to compute on modern computers, but it will also still share many of its characteristics with p
→
linear
(
t
) and allow trajectory planning to be done with relative ease. The remainder of this paper mainly revolves around exploring these and other properties of p
→
(
t
) together with their implications.


### 3.3. A Qualitative Comparison

Asp
→
(
t
) is only an approximation, it will differ from the linear drag model’s trajectory function it is based on, as well as from the results of other models. For comparison purposes, the trajectories that follow from launching three different sport balls using four different models are plotted side by side in Figure 3. All balls are launched at an 45° angle at 50 m/s on a windless earthly day (g
=
9
.
81 m/s2). The different drag models are calibrated to respect the respective ball’s terminal velocities (except for the dragless model, which always has an infinite terminal velocity).Figure 3
The trajectories resulting from four models for three different ball types.Each of the alternating thick and thin segments in the trajectories shown in Figure3 (and in all other trajectory plots in this paper, for that matter) represents a projectile’s movement over a period of exactly one second, making it possible to not only compare the shapes of the trajectories but also their local speeds. The results from the novel p
→
(
t
) function are plotted in black, the results from the linear drag model p
→
linear
(
t
) are plotted in blue, and the results of the physically (most) correct quadratic drag model simulation p
→
quadratic
(
t
) are plotted in green. Lastly, the red parabola p
→
dragless
(
t
) represents the trajectory of each of the three balls in a perfect vacuum (in other words, when there is no drag and v
→
∞’s length goes to infinite).When comparing the trajectories forp
→
(
t
) to the results of the two more accurate drag models, they are certainly different but they still reasonably mimic these in look, feel, and properties. Consequently, the proposed model is physically at least quite plausible and is probably accurate enough for most computer game purposes. Furthermore, in some cases, p
→
(
t
) is actually closer to p
→
quadratic
(
t
) than p
→
linear
(
t
), making it in these cases arguably even physically more accurate than the model it is approximating to. Lastly, the trajectories for all three drag models perfectly approach p
→
dragless
(
t
) when v
→
∞’s length goes to infinite.


### 3.4. Exploring Some ofp
(
t
)’s Properties

The functionp
→
(
t
) given by (4) can be factored, solved, and parameterized in many different ways. For example, basic algebra allows it to be written as p
→
(
t
)
=
a
v
→
0
+
b
→ as well, where a
=
t
/
(
1
+
k
t
) and b
→
=
v
→
∞
(
t
-
a
)
+
p
→
0. Note that in this form, the initial velocity is separated from all other factors, and it becomes immediately clear that p
→
(
t
) is a linear function in terms of v
→
0. This implies that when launching multiple objects at some t
=
0 with all properties equal except for the initial velocity, each of these objects has the same value for a and for b
→. This feature may be exploited in particle explosion systems on modern GPUs, for example, requiring only one evaluation of a and b
→ per frame per explosion (layer) on the CPU and one MAD (multiply-and-add) GPU instruction per particle per frame to calculate each particle’s position.The linearity ofp
→
(
t
) in terms of the initial velocity v
→
0 can also be used for many other purposes. For example, in Figure 4, p
→
(
t
) is used to calculate the green and blue positions for two different “extreme” initial velocities, which are interpreted as the top-left and bottom-right positions of a textured rectangle orquad. Note that the (signed) size of the quad is thus simply p
→
green
(
t
)
-
p
→
blue
(
t
)
=
(
a
v
→
0
,
green
+
b
→
)
-
(
a
v
→
0
,
blue
+
b
→
)
=
a
(
v
→
0
,
green
-
v
→
0
,
blue
). Furthermore, all the bilinearly interpolated texels within this quad, including the red one, will move over p
→
(
t
) trajectories themselves as well by virtue of the linearity in terms of the initial velocities. In other words, each texel will follow a p
→
(
t
) trajectory with some initial velocity that is interpolated bilinearly between the different extreme initial velocities, as if the individual texels themselves are under direct control of a physics simulation. Consequently, it should be physically plausible to use the above to scale sprites and billboards of, for example, simple smoke (for which v
→
∞ would typically be upwards to simulate positive buoyancy), explosion debris, and fireworks.Figure 4
The trajectory of a textured quad.Many other useful properties can easily be derived fromp
→
(
t
) as well. For example, the velocity v
→
(
t
) of p
→
(
t
) is
(5)
v
→
(
t
)
=
∂
p
→
(
t
)
∂
t
=
v
→
0
+
k
t
(
2
+
k
t
)
v
→
∞
(
1
+
k
t
)
2
.The nonnegative timet
top
,
n
→ at which the trajectory hits its maximum in the direction of a given unit-length vector n
→ can be found by solving v
→
(
t
)
·
n
→
=
0 for t assuming that the direction to find the maximum in is pointing away from v
→
∞ (i.e., v
→
∞
·
n
→
<
0). If that assumption is false, then the top will be at time t
=
0. The solution to both cases is summarized by the following formula:
(6)
t
top
,
n
→
=
1
-
min
⁡
(
0
,
(
v
→
0
·
n
→
)
/
(
v
→
∞
·
n
→
)
)
-
1
k
.Thist
top
,
n
→ may be used with (4) to get the trajectory’s maximum position in the n
→ direction. Note that when n
→ is axis-aligned, the dot products in (6) can be optimized away. For example, when n
→ is equal to the +y axis, then t
top
,
n
→ becomes
(7)
t
top
,
y
=
1
-
min
⁡
(
0
,
v
→
0
,
y
/
v
→
∞
,
y
)
-
1
k
.See Figure5 for an example of v
→
(
t
) and p
→
(
t
top
,
y
).Figure 5
The velocity andy-top of a trajectory.


## 3.1. Ballistic Parameters

Before presenting the proposed trajectory function and its properties, the necessary ballistic parameters will be defined here first. Starting with a general note, vector variables in this paper are always distinguished from scalar variables by the→ symbol above their names. Also, the length of any vector a
→ is denoted as ∥
a
→
∥, which is equal to a
→
·
a
→, where · is the dot product.A ballistic object is assumed to be launched from the initial positionp
→
0 with the initial velocity v
→
0 at time t
=
0. Furthermore, the object will be travelling through a medium (e.g., the air) which itself travels at the fixed (wind) velocity v
→
medium. It is also pulled by gravity at the fixed gravitational acceleration g, which is roughly 9.81 m/s2 for “earthly” applications. The amount of drag k while moving through the medium is defined as follows.
(1)
k
=
1
2
g
∥
v
→
terminal
∥
.Here,v
→
terminal is the invariant terminal velocity relative to v
→
medium that is reached eventually as the forces of gravity and air resistance finally cancel each other out. On earth, that is equivalent to saying that v
→
terminal is the fixed velocity that is approached when the object is dropped from an enormous height on a windless day. Lastly, the absolute terminal velocity v
→
∞, being the absolute velocity approached when time t goes to ∞, is therefore
(2)
v
→
∞
=
v
→
terminal
+
v
→
medium
.Together,p
→
0, v
→
0, v
→
∞, and k uniquely define a trajectory in the proposed model. To give a real-world example of the parameters defined above, suppose a tennis ball with a terminal velocity of 30 m/s is served at 50 m/s at a 30° angle from the 381 m high roof of the Empire State building into a 10 m/s horizontal wind in 2D. Then, p
→
0
=
[
0
,
381
], v
→
0
=
[
50
cos
⁡
(
3
0
°
)
,
50
sin
(
3
0
°
)
]
≈
[
43.30
,
25
], k
=
0.1635, and v
→
∞
=
[
-
10
,
-
30
]. The trajectory resulting from these values is shown in Figure 1.Figure 1
Serving a tennis ball from the Empire State Building upwind.


## 3.2. Deriving the Trajectory Function

In terms of the parameters defined above, the differential equation for the exact linear drag model has the following analytic solution:(3)
p
→
linear
(
t
)
=
1
2
k
(
v
→
0
-
v
→
∞
)
(
1
-
e
-
2
k
t
)
+
v
→
∞
t
+
p
→
0
.This function calculates the 2D or 3D positionp
→
linear on a trajectory at time t where t
≥
0. The above function is far from new and will not be explained here in detail, as it has already been covered in many textbooks [4], occasionally even targeting game developers in particular [9] (albeit with slightly different notation and parameter definitions).The function above will not be used directly in this paper. Instead, it will be approximated by substituting its exponential functione
x with the first degree rational function (
2
+
x
)
/
(
2
-
x
) shown in Figure 2. One of the reasons for selecting this approximation over all possible other approximations to e
x is that it has a value, first derivative and second derivative that match those of e
x at x
=
0. This means that it approximates e
x near x
=
0 well and therefore will guarantee a good approximation of p
→
linear
(
t
) near t
=
0. Furthermore, the first derivative of (
2
+
x
)
/
(
2
-
x
) monotonically decreases from 1 to 0 as x tends from 0 to -
∞, similar to the first derivative of e
x itself. When used to approximate e
x in p
→
linear
(
t
), this property will cause the initial velocity to monotonically converge to the terminal velocity over time. Note that no polynomial approximation of e
x has this specific property, and of all possible rational functions that do possess the above properties, the proposed approximation is the simplest and thus the most efficient. Lastly, its inverse is also a first degree rational function, resulting in relatively simple algebraic solutions for all (otherwise algebraic) equations that use it to approximate e
x.Figure 2
Approximatinge
x using a simple rational function.When the exponential function inp
→
linear
(
t
) is substituted by the rational approximation, the following function is the result
(4)
p
→
(
t
)
=
(
v
→
0
+
k
t
v
→
∞
)
t
1
+
k
t
+
p
→
0
.Because of the aforementioned properties,p
→
(
t
) will not only be more efficient to compute on modern computers, but it will also still share many of its characteristics with p
→
linear
(
t
) and allow trajectory planning to be done with relative ease. The remainder of this paper mainly revolves around exploring these and other properties of p
→
(
t
) together with their implications.


## 3.3. A Qualitative Comparison

Asp
→
(
t
) is only an approximation, it will differ from the linear drag model’s trajectory function it is based on, as well as from the results of other models. For comparison purposes, the trajectories that follow from launching three different sport balls using four different models are plotted side by side in Figure 3. All balls are launched at an 45° angle at 50 m/s on a windless earthly day (g
=
9
.
81 m/s2). The different drag models are calibrated to respect the respective ball’s terminal velocities (except for the dragless model, which always has an infinite terminal velocity).Figure 3
The trajectories resulting from four models for three different ball types.Each of the alternating thick and thin segments in the trajectories shown in Figure3 (and in all other trajectory plots in this paper, for that matter) represents a projectile’s movement over a period of exactly one second, making it possible to not only compare the shapes of the trajectories but also their local speeds. The results from the novel p
→
(
t
) function are plotted in black, the results from the linear drag model p
→
linear
(
t
) are plotted in blue, and the results of the physically (most) correct quadratic drag model simulation p
→
quadratic
(
t
) are plotted in green. Lastly, the red parabola p
→
dragless
(
t
) represents the trajectory of each of the three balls in a perfect vacuum (in other words, when there is no drag and v
→
∞’s length goes to infinite).When comparing the trajectories forp
→
(
t
) to the results of the two more accurate drag models, they are certainly different but they still reasonably mimic these in look, feel, and properties. Consequently, the proposed model is physically at least quite plausible and is probably accurate enough for most computer game purposes. Furthermore, in some cases, p
→
(
t
) is actually closer to p
→
quadratic
(
t
) than p
→
linear
(
t
), making it in these cases arguably even physically more accurate than the model it is approximating to. Lastly, the trajectories for all three drag models perfectly approach p
→
dragless
(
t
) when v
→
∞’s length goes to infinite.


## 3.4. Exploring Some ofp
(
t
)’s Properties

The functionp
→
(
t
) given by (4) can be factored, solved, and parameterized in many different ways. For example, basic algebra allows it to be written as p
→
(
t
)
=
a
v
→
0
+
b
→ as well, where a
=
t
/
(
1
+
k
t
) and b
→
=
v
→
∞
(
t
-
a
)
+
p
→
0. Note that in this form, the initial velocity is separated from all other factors, and it becomes immediately clear that p
→
(
t
) is a linear function in terms of v
→
0. This implies that when launching multiple objects at some t
=
0 with all properties equal except for the initial velocity, each of these objects has the same value for a and for b
→. This feature may be exploited in particle explosion systems on modern GPUs, for example, requiring only one evaluation of a and b
→ per frame per explosion (layer) on the CPU and one MAD (multiply-and-add) GPU instruction per particle per frame to calculate each particle’s position.The linearity ofp
→
(
t
) in terms of the initial velocity v
→
0 can also be used for many other purposes. For example, in Figure 4, p
→
(
t
) is used to calculate the green and blue positions for two different “extreme” initial velocities, which are interpreted as the top-left and bottom-right positions of a textured rectangle orquad. Note that the (signed) size of the quad is thus simply p
→
green
(
t
)
-
p
→
blue
(
t
)
=
(
a
v
→
0
,
green
+
b
→
)
-
(
a
v
→
0
,
blue
+
b
→
)
=
a
(
v
→
0
,
green
-
v
→
0
,
blue
). Furthermore, all the bilinearly interpolated texels within this quad, including the red one, will move over p
→
(
t
) trajectories themselves as well by virtue of the linearity in terms of the initial velocities. In other words, each texel will follow a p
→
(
t
) trajectory with some initial velocity that is interpolated bilinearly between the different extreme initial velocities, as if the individual texels themselves are under direct control of a physics simulation. Consequently, it should be physically plausible to use the above to scale sprites and billboards of, for example, simple smoke (for which v
→
∞ would typically be upwards to simulate positive buoyancy), explosion debris, and fireworks.Figure 4
The trajectory of a textured quad.Many other useful properties can easily be derived fromp
→
(
t
) as well. For example, the velocity v
→
(
t
) of p
→
(
t
) is
(5)
v
→
(
t
)
=
∂
p
→
(
t
)
∂
t
=
v
→
0
+
k
t
(
2
+
k
t
)
v
→
∞
(
1
+
k
t
)
2
.The nonnegative timet
top
,
n
→ at which the trajectory hits its maximum in the direction of a given unit-length vector n
→ can be found by solving v
→
(
t
)
·
n
→
=
0 for t assuming that the direction to find the maximum in is pointing away from v
→
∞ (i.e., v
→
∞
·
n
→
<
0). If that assumption is false, then the top will be at time t
=
0. The solution to both cases is summarized by the following formula:
(6)
t
top
,
n
→
=
1
-
min
⁡
(
0
,
(
v
→
0
·
n
→
)
/
(
v
→
∞
·
n
→
)
)
-
1
k
.Thist
top
,
n
→ may be used with (4) to get the trajectory’s maximum position in the n
→ direction. Note that when n
→ is axis-aligned, the dot products in (6) can be optimized away. For example, when n
→ is equal to the +y axis, then t
top
,
n
→ becomes
(7)
t
top
,
y
=
1
-
min
⁡
(
0
,
v
→
0
,
y
/
v
→
∞
,
y
)
-
1
k
.See Figure5 for an example of v
→
(
t
) and p
→
(
t
top
,
y
).Figure 5
The velocity andy-top of a trajectory.


## 4. Planning in World Space: Hitting a Target at a Given Time

When a projectile needs to hit some given target position,p
→
(
t
) can be used to solve or “plan” the initial velocity that leads to precisely hitting this target at some given future point in time. To be more specific, when trying to hit some position r
→ at the given time t
r, the solution is found by solving p
→
(
t
r
)
=
r
→ for v
→
0, which results in the following formula:
(8)
v
→
0
=
k
(
r
→
-
p
→
0
-
v
→
∞
t
r
)
+
r
→
-
p
→
0
t
r
.In Figure6, this formula is used to plan the trajectories to six different target r
→ positions, all taking exactly ten seconds to reach their target. In other words, t
r
=
10, which results in having exactly five thick and five thin segments per trajectory in this figure.Figure 6
Planning to hit all six different target positions in exactly ten seconds.One interesting property of this function is that thex and y components (or the x, y, and z components in the 3D case) of v
→
0 are completely independent from each other. As a direct consequence, similar projectiles that target the same horizontal distance but a different height will always move at the same horizontal speed over the whole trajectory and vice versa. Also, trajectories for targets at equal height but different horizontal distances will all have the same top height, as can be observed in Figure 6 as well.Note that the given time parametert
r could be any positive value, including the one that is dependent on another function of r
→. One simple example of such a function is t
r
=
∥
r
→
-
p
→
0
∥
/
v
straight, where v
straight is a given average speed over the straight line from p
→
0 to r
→. In Figure 7, this is illustrated using v
straight
=
25 m/s.Figure 7
Planning trajectories to six targets, taking the exact same time as travelling with 25 m/s over the straight line from initial to target position.A function fort
r could be made arbitrarily complex. The time planners that will be explored in Section 6 are examples of moderately complex functions, solving different additional constraints for t
r given r
→. But even more complex planners would be necessary if t
r and r
→ were to be dependent on each other. This can happen, for example, when a moving target needs to be hit while planning to launch at a fixed speed, making t
r influence the prediction of the future target position r
→, which influences t
r again. These relationships are not explored in detail in this paper, but it is worth mentioning that some of these problems may be solved iteratively by starting with a rough estimate for t
r and then letting it converge to the right solution by repeatedly going from t
r to r
→ and from r
→ to an improved t
r. These iterations could possibly even be spread over multiple frames to amortize costs, for example, improving accuracy with each new frame.


## 5. The Principal Frame of Reference and Its Properties

Most planners forp
→
(
t
) are still to be presented. However, as these planners depend on a special frame of reference to keep the required planner math as simple as possible, this frame of reference will be covered here first.Inspectingp
→
(
t
) reveals that the function always outputs p
→
0 plus a linear combination (i.e., a weighted sum) of v
→
0 and v
→
∞. Geometrically, this implies that all trajectories, even with wind coming from any 3D direction, are guaranteed to lie on a plane spanned by v
→
0 and v
→
∞ which passes through p
→
0. Consequently, p
→
(
t
) may be rewritten as
(9)
p
→
(
t
)
=
x
(
t
)
X
→
+
y
(
t
)
Y
→
+
p
→
0
,
where X
→ and Y
→, respectively, describe an orthonormal tangent and bitangent direction inworld space of the plane over which the projectile moves. Within this plane’s 2D frame of reference, the vector [
x
(
t
)
,
y
(
t
)
] describes the movement on the trajectory over time relative to p
→
0 and in terms of this alternative x axis (i.e., the tangent) and y axis (i.e., the bitangent).Put into terms perhaps more familiar to computer graphics and game developers, the 2D function[
x
(
t
)
,
y
(
t
)
] is like a procedural UV coordinate that defines a projectile’s location on an unwrapped plane which maps UV [
0
,
0
] to p
→
0 and has an orthonormal tangent vector X
→ and bitangent vector Y
→. All this is shown in Figure 8, where the trajectory p
→
(
t
) is visualized as the intersection between the described plane (on which is thus lies) and another curved surface.Figure 8
A 3D trajectory with wind at an angle, and its decomposition on the principal plane (in blue).Although there is an infinite amount of ways to define theX
→ and Y
→ vectors, the following definitions are used in this paper for their particularly useful properties:
(10)
X
→
=
d
→
-
(
d
→
·
Y
→
)
Y
→
∥
d
→
-
(
d
→
·
Y
→
)
Y
→
∥
,
Y
→
=
-
v
→
∞
v
∞
.Here,v
∞ is used as a shorthand for ∥
v
→
∞
∥. And d
→
=
v
→
0 if v
→
0 is already known. In the case that v
→
0 is not (yet) known, any position relative to p
→
0 known to be lying on the trajectory may be used for d
→ instead. For example, use d
→
=
r
→
-
p
→
0 when targeting the position r
→. Note that in the case that v
→
∞ and d
→ are collinear, the trajectory can be described by movement solely in the Y
→ direction. To still get a valid 2D basis in that case, an arbitrary vector that is noncollinear should be used for d
→ instead.The frame of reference defined byp
→
0, X
→, and Y
→ is what will be called the trajectory’sprincipal space. It is called that because this space allows the math describing the trajectory to be decomposed into a particularly compact and well-behaved form. In particular, wind and gravity do not affect movement over the principal X
→ axis at all, but solely over the principal Y
→ axis. This allows the function x
(
t
) in (9) to be simplified as will be shown and put to good use soon.The second useful property of this principal space is that it guarantees thatv
x
≥
0. That is, all initial and in-flight velocities expressed in principal space are guaranteed to have nonnegative values on the x axis, causing ballistic objects to never move to or be on the negative x side of this space, even though they obviously can still move in any direction in world space. Similarly, targets are never on the negative x side in principal space either. The advantage of this property is that it once again will allow for further simplifications in some of the planner math that is still to be discussed.The third useful property of this particular space is that it is easy to convert from world space to principal space and back asX
→ and Y
→ are orthonormal. Converting from any world space position p
→ to the position [
x
,
y
] in principal space and vice versa can simply be done using (11) and (12), respectively as follows:
(11)
[
x
,
y
]
=
[
(
p
→
-
p
→
0
)
·
X
→
,
(
p
→
-
p
→
0
)
·
Y
→
]
,
(12)
p
→
=
x
X
→
+
y
Y
→
+
p
→
0
.Similarly, it is possible to efficiently convert from any world space velocityv
→ to the principal velocity [
v
x
,
v
y
] and back using (13) and (14), respectively as follows:
(13)
[
v
x
,
v
y
]
=
[
v
→
·
X
→
,
v
→
·
Y
→
]
,
(14)
v
→
=
v
x
X
→
+
v
y
Y
→
.Starting a new notational convention here for clarity, vector names (i.e., variables decorated with a→ symbol) are only used for variables in world space, while variables in principal space never use this decoration and always represent individual scalar quantities. So, for example, v
→
x is the x component of the vector representing the world space velocity v
→, while v
x is a scalar representing a velocity in the x direction in principal space.Now that the frame of reference itself has been covered, it is possible to define the two scalar functions that make up the principal space trajectory function[
x
(
t
)
,
y
(
t
)
]:
(15)
x
(
t
)
=
v
0
,
x
t
1
+
k
t
,
(16)
y
(
t
)
=
(
v
0
,
y
-
k
t
v
∞
)
t
1
+
k
t
.These functions are derived by transformingp
→
(
t
) into this space using (11). Note that gravity and wind do indeed not affect movement over the x direction in this space. And [
x
(
0
)
,
y
(
0
)
] is equal to [
0
,
0
], which means that trajectories in principal space always start at the origin (while starting at p
→
0 in world space).The simpler formula forx
(
t
) makes it possible to uniquely invert the function to get the time t at which the x component of a certain position will be reached given the horizontal initial velocity. The solution is as follows:
(17)
t
=
x
v
0
,
x
-
k
x
.Here,0
≤
x
≤
v
0
,
x
/
k, as that is the valid range of x
(
t
) for t
≥
0 as defined by (15). By plugging (17) into (16), the following explicit y-for-x relationship is found:
(18)
y
(
x
)
=
(
k
v
∞
x
/
(
k
x
-
v
0
,
x
)
+
v
0
,
y
)
x
v
0
,
x
.This function always has exactly oney value for each valid x value, which would not necessarily be true for an explicit trajectory function in any other space. This property is demonstrated in Figure 9, showing a trajectory that is equivalent to the trajectory shown in Figure 1 but which is now plotted using y
(
x
) in principal space, perfectly overlaying the original trajectory when mapped back into world space.Figure 9
The same Empire State Building trajectory example as from Figure1 but seen from the principal frame of reference.For completeness, the principal space counterparts of the world space properties described by (5), (6), and (7) are given here as well. Consider
(19)
v
x
(
t
)
=
v
0
,
x
(
1
+
k
t
)
2
,
v
y
(
t
)
=
v
0
,
y
-
k
t
(
2
+
k
t
)
v
∞
(
1
+
k
t
)
2
,
t
top
,
n
→
=
1
+
max
⁡
(
0
,
[
v
0
,
x
,
v
0
,
y
]
·
[
n
x
,
n
y
]
/
(
v
∞
n
y
)
)
-
1
k
t
top
,
y
=
1
+
max
⁡
(
0
,
v
0
,
y
/
v
∞
,
y
)
-
1
k
.Lastly, the local slope in principal space in terms of time (i.e.,v
y
(
t
)
/
v
x
(
t
)) and in terms of x (i.e., ∂
y
(
x
)
/
∂
x) is as follows:
(20)
v
y
(
t
)
v
x
(
t
)
=
v
0
,
y
-
k
t
(
2
+
k
t
)
v
∞
v
0
,
x
,
∂
y
(
x
)
∂
x
=
v
0
,
y
+
v
∞
v
0
,
x
-
v
0
,
x
v
∞
(
k
x
-
v
0
,
x
)
2
.


## 6. Planning in Principal Space

As the planners in this section all depend on the properties of trajectories in principal space, the most relevant properties are briefly repeated here. Per definition, any ballistic object in principal space is launched from the origin, any target has a nonnegativex component, and the combined effect of gravity and wind results in a v
∞ value that is exactly in the -
y direction. As the planners will expect their parameters to be specified in principal space, the parameters of any world space problem need to be converted to this space before they can be used. To recap the necessary steps (assuming the problem involves hitting some target position r
→), start by defining the actual principal space’s X
→ and Y
→ axes using (10) with d
→
=
r
→
-
p
→
0. Next, convert r
→ (or any other requested position) to principal space using (11) to get [
r
x
,
r
y
].All planners covered here will return the exact timet
r at which the target [
r
x
,
r
y
] must be hit to meet the planner’s given constraints. To get the actual initial velocity in principal space that leads to hitting [
r
x
,
r
y
] at this t
r, both p
x
(
t
r
)
=
r
x and p
y
(
t
r
)
=
r
y need to be solved for [
v
0
,
x
,
v
0
,
y
], which can be done using the following two formulas:
(21)
v
0
,
x
=
r
x
(
k
+
1
t
r
)
,
(22)
v
0
,
y
=
r
y
(
k
+
1
t
r
)
+
k
v
∞
t
r
.To get the initial velocity in world space from these, it is possible to convert[
v
0
,
x
,
v
0
,
y
] to v
→
0 using (14). But v
→
0 may also be directly calculated from t
r and r
→ through (8). Now that it is clear how to make use of principal space planners in general, the actual planners are presented.


### 6.1. Hitting the Target Given Another Position to Pass Through

A trajectory can be planned to pass through both position[
q
x
,
q
y
] and through target [
r
x
,
r
y
] by solving y
(
q
x
)
=
q
y and y
(
r
x
)
=
r
y for v
0
,
x, and using (17) on r
x and v
0
,
x to get t
r. This specific form of planning may be useful to shoot through a hole or exactly over an object at [
q
x
,
q
y
] to hit [
r
x
,
r
y
], example. The solution to at what time the position [
r
x
,
r
y
] needs to be hit is as follows:
(23)
t
r
=
s
+
s
(
s
+
2
r
x
k
(
r
x
-
q
x
)
)
,
where s
=
(
r
x
q
y
-
r
y
q
x
)
/
2
v
∞
q
x. Note that the line from the origin to the target position with the smallest x must be at least as steep as the line from the origin to the position with the largest x for a valid trajectory (and thus a real solution) to exist. That is, q
y
/
q
x
≥
r
y
/
r
x if q
x
<
r
x, and r
y
/
r
x
≥
q
y
/
q
x if q
x
>
r
x.  In Figure 10, two of the six targets do not meet this requirement, explaining why there are only four trajectories there.Figure 10
Planning trajectories throughQ and six different R targets.


### 6.2. Hitting the Target While Touching a Line

When looking for the timet
r at which a trajectory passes through [
r
x
,
r
y
] and touches the line y
=
a
x
+
b in principal space, ∂
y
(
x
)
/
∂
x
=
a has to be solved for x first, which can then be used to solve y
(
x
)
=
a
x
+
b for v
0
,
x. Both v
0
,
x and r
x can then be used again with (17) to get t
r. The solution may be written as follows.                              (24)
t
r
=
s
+
s
2
+
k
v
∞
(
a
r
x
-
r
y
)
k
v
∞
,
where s
=
b
k
v
∞
+
(
1
/
2
)
k
(
a
r
x
-
r
y
+
b
). Note that a real solution can only exist when b
≥
0 and a
r
x
+
b
≥
r
y. That is, the line must always pass through or be above the initial and the target position.In the example presented in Figure11, the line to be touched is chosen to be horizontal, leading to a specification of the trajectories’ vertical tops in principal space. But it is also possible to find the t
r that leads to hitting a top defined in another space. For example, to let a trajectory’s top touch the world-space plane with normal n
→ and through point q
→, y
=
a
x
+
b becomes the line in principal space that describes the intersection between this plane and the trajectory’s principal plane. In that case, a
=
-
(
X
→
·
n
→
)
/
(
Y
→
·
n
→
) and  b
=
(
(
q
→
-
p
→
0
)
·
n
→
)
/
(
Y
→
·
n
→
). As always, when n
→ is axis-aligned, the dot products can be optimized away. For example, when only interested in trajectories exactly hitting a world space height h at their tops (in the +y direction), this simplifies to a
=
-
X
→
y
/
Y
→
y and b
=
(
h
-
p
→
0
,
y
)
/
Y
→
y.Figure 11
Planning trajectories through six different targets, touching a line at each top.The principal space slopea can also be calculated from the slope w of a world-space elevation angle θ
world (i.e., w
=
tan
(
θ
world
)) by using the conversion formula a
=
(
sign
⁡
(
Y
→
y
)
w
s
(
X
→
y
2
+
Y
→
y
2
)
-
w
2
-
s
X
→
y
Y
→
y
)
/
(
s
Y
→
y
2
-
w
2
), where s
=
1
+
w
2 and sign
⁡
(
x
)
=
[
x
≥
0
]
-
[
x
<
0
]. This a is only valid if an equivalent elevation with a positive x component in principal space exists, which is the case if s
(
X
→
y
2
+
Y
→
y
2
)
-
w
2
≥
0 and [
1
,
w
]
·
[
1
-
X
→
y
2
,
X
→
y
]
>
0. This world-space elevation conversion may be particularly useful when used together with the next two principal-space planners.


### 6.3. Hitting the Target Given the Initial Slope

This subsection is about finding the timet
r at which a projectile will hit position [
r
x
,
r
y
] while being launched at slope a. Planning this way may be useful when there is control over the projectile’s initial speed but not over its direction (e.g., for some weapon mounted on a fixed rig).This problem is actually simply a special case of the previous planner, whereb
=
0, meaning that the problem is equivalent to finding the trajectory that touches the line y
=
a
x. After substituting b in (24) and applying some basic algebra, the solution may be written more compactly as follows:
(25)
t
r
=
s
+
s
(
s
+
2
v
∞
)
k
v
∞
,
where s
=
(
1
/
2
)
k
(
a
r
x
-
r
y
). Obviously, the slope a has to be steeper than the line from the origin to the target position (i.e., a
>
r
y
/
r
x) for a solution to exist, which is the case for all but one target position in Figure 12.Figure 12
Planning trajectories through six different targets, all starting with the slope 3/4.


### 6.4. Hitting the Target Given the Target Slope

Similarly, it is possible to hit[
r
x
,
r
y
] given the exact slope a at the target position. This type of planning allows for exact control over the angle at which a target is hit. Again, this can be seen as a special case of the planner from Section 6.2, using r
y
-
a
r
x as the value for b. When substituted and simplified, this results in the following more direct formula:
(26)
t
r
=
r
y
-
a
r
x
k
v
∞
.For this particular planner, the slopea needs to be less than the slope of the line from the origin to the target position (i.e., a
<
r
y
/
r
x). Consequently, only five of the six targets have a valid trajectory in Figure 13.Figure 13
Planning trajectories through six different targets, all ending at the slope −3/4.


### 6.5. Hitting the Target Given the Arc Height or “Curviness”

The timet
r can also be calculated for a target position [
r
x
,
r
y
] and an “arc height” b. Here, b is defined as the maximum difference in the y direction between the trajectory and the straight line from origin to target. Equivalently, b may be interpreted as the height in principal space of the smallest parallelogram containing the whole trajectory, as shown in Figure 14. This problem is another special case of the “line touching” problem and can be solved by substituting r
y
/
r
x for a in (24). After applying some basic algebra, the resulting formula may be written as follows:
(27)
t
r
=
b
v
∞
+
2
b
k
v
∞
.Figure 14
Planning trajectories through six differentRs, all with an arc height of 50.This function is particularly intuitive to plan with when(28)
b
=
h
r
x
2
+
r
y
2
,
where h defines the “curviness” of the trajectory. For example, h
=
0.01 always leads to a low arc for any target position, while h
=
0.5 always leads to a fairly high arc.


### 6.6. Hitting the Target with (Almost) Minimal Effort

A trajectory can also be planned to hit a target at[
r
x
,
r
y
] with the smallest initial speed (and thus the least amount of energy) possible. Planning a trajectory this way requires finding the positive time-to-target t
r which solves ∂
(
v
0
,
x
2
+
v
0
,
y
2
)
/
∂
t
r
=
0, where v
0
,
x and v
0
,
y are defined by (21) and (22), respectively. This equation can be expanded into the following form:
(29)
(
k
2
v
∞
2
)
t
r
4
+
(
k
2
v
∞
r
y
)
t
r
3
-
k
(
r
x
2
+
r
y
2
)
t
r
-
(
r
x
2
+
r
y
2
)
=
0
.Like all quartic equations, solving this in a closed form is possible but difficult to do robustly [10]. In practice, quartic equations are typically solved for any or all of their roots by generic iterative root solvers [11]. But by exploiting domain-specific knowledge, it is also possible to implement a specialized iterative solver for (29) that is guaranteed to efficiently converge to the right root directly. One possible implementation is presented as a C++ function calledGetTimeToTargetRWithMinimalInitialSpeed() in Algorithm 1. There, the equation is first transformed into an equivalent but more well-behaved strictly convex quartic function for which a conservative initial guess for t (or rather, u) is calculated, which is then refined using multiple (but typically less than a total of six) modified or normal conservative Newton’s method iterations. This implementation has been carefully crafted with both robustness and efficiency in mind. As with any implementation, numeric precision can become an issue when using extreme values, but results for practical ranges are typically within a fewfloat epsilons of the exact value. See the comments in the implementation itself for more details.Alternatively, when only a rough approximation of the minimal effort solution is needed, the simpler “curviness” planner from (27) and (28) could be used with h
=
1
/
4. This approximation is fairly accurate for larger values of v
∞, as it actually converges perfectly to the exact solution when v
∞ goes to infinity. But for high friction scenarios, the difference between the exact method and the approximation becomes quite noticeable. To get an idea of the size of the error for a medium friction scenario, the trajectories resulting from the exact method and the approximation are shown side by side in Figure 15 for the case of g
=
9.81m/s2 and v
∞
=
25m/s.Figure 15
Planning trajectories through six differentRs while minimizing the initial velocity using the exact method (in black) and the approximate h
=
1
/
4 method (in red) in a medium-fiction scenario.


### 6.7. Hitting a Target Given the Initial Speed

The last planning algorithm covered in this paper solves for the timet
r required to hit the target [
r
x
,
r
y
] given the projectile’s exact initial speed s. The solution is found by solving v
0
,
x
2
+
v
0
,
y
2
=
s
2 for the positive time t
r to position [
r
x
,
r
y
], where v
0
,
x and v
0
,
y are once again defined by (21) and (22), respectively. This particular equation may be expanded into the following quartic function:
(30)
(
k
2
v
∞
2
)
t
r
4
+
(
2
k
2
v
∞
r
y
)
t
r
3
+
(
k
2
r
x
2
+
k
2
r
y
2
+
2
k
v
∞
r
y
-
s
2
)
t
r
2
+
2
k
(
r
x
2
+
r
y
2
)
t
r
+
(
r
x
2
+
r
y
2
)
=
0
.Note that whens is smaller than the minimal initial speed to hit [
r
x
,
r
y
] (i.e., the sought solution of (29)), the problem has no valid solutions. But when s is larger than that, it will always have exactly two valid solutions. In that case, the smaller of the two t
r values represents the time to hit the target with a low arc, while the larger t
r root is the solution for a high arc. See Figure 16 for an example.Figure 16
Planning trajectories through six differentRs, all starting at 100 m/s. Each of the five reachable targets has a trajectory with a (red) low arc and a (black) high arc.Using a similar approach as in the previous section, a specialized iterative solver can be implemented to solve this particular quartic function. The C++ functionGetTimeToTargetRGivenInitialSpeedS() in Algorithm 2 solves the equation for either the high or low arc root and returns the resulting t
r or returns 0 if no solution exists. This implementation has also been written with robustness, efficiency, and accuracy for a wide range of parameters in mind. Tests showed that the procedure typically requires about six (modified) Newton’s method iterations in total to converge to almost fullfloat precision.


## 6.1. Hitting the Target Given Another Position to Pass Through

A trajectory can be planned to pass through both position[
q
x
,
q
y
] and through target [
r
x
,
r
y
] by solving y
(
q
x
)
=
q
y and y
(
r
x
)
=
r
y for v
0
,
x, and using (17) on r
x and v
0
,
x to get t
r. This specific form of planning may be useful to shoot through a hole or exactly over an object at [
q
x
,
q
y
] to hit [
r
x
,
r
y
], example. The solution to at what time the position [
r
x
,
r
y
] needs to be hit is as follows:
(23)
t
r
=
s
+
s
(
s
+
2
r
x
k
(
r
x
-
q
x
)
)
,
where s
=
(
r
x
q
y
-
r
y
q
x
)
/
2
v
∞
q
x. Note that the line from the origin to the target position with the smallest x must be at least as steep as the line from the origin to the position with the largest x for a valid trajectory (and thus a real solution) to exist. That is, q
y
/
q
x
≥
r
y
/
r
x if q
x
<
r
x, and r
y
/
r
x
≥
q
y
/
q
x if q
x
>
r
x.  In Figure 10, two of the six targets do not meet this requirement, explaining why there are only four trajectories there.Figure 10
Planning trajectories throughQ and six different R targets.


## 6.2. Hitting the Target While Touching a Line

When looking for the timet
r at which a trajectory passes through [
r
x
,
r
y
] and touches the line y
=
a
x
+
b in principal space, ∂
y
(
x
)
/
∂
x
=
a has to be solved for x first, which can then be used to solve y
(
x
)
=
a
x
+
b for v
0
,
x. Both v
0
,
x and r
x can then be used again with (17) to get t
r. The solution may be written as follows.                              (24)
t
r
=
s
+
s
2
+
k
v
∞
(
a
r
x
-
r
y
)
k
v
∞
,
where s
=
b
k
v
∞
+
(
1
/
2
)
k
(
a
r
x
-
r
y
+
b
). Note that a real solution can only exist when b
≥
0 and a
r
x
+
b
≥
r
y. That is, the line must always pass through or be above the initial and the target position.In the example presented in Figure11, the line to be touched is chosen to be horizontal, leading to a specification of the trajectories’ vertical tops in principal space. But it is also possible to find the t
r that leads to hitting a top defined in another space. For example, to let a trajectory’s top touch the world-space plane with normal n
→ and through point q
→, y
=
a
x
+
b becomes the line in principal space that describes the intersection between this plane and the trajectory’s principal plane. In that case, a
=
-
(
X
→
·
n
→
)
/
(
Y
→
·
n
→
) and  b
=
(
(
q
→
-
p
→
0
)
·
n
→
)
/
(
Y
→
·
n
→
). As always, when n
→ is axis-aligned, the dot products can be optimized away. For example, when only interested in trajectories exactly hitting a world space height h at their tops (in the +y direction), this simplifies to a
=
-
X
→
y
/
Y
→
y and b
=
(
h
-
p
→
0
,
y
)
/
Y
→
y.Figure 11
Planning trajectories through six different targets, touching a line at each top.The principal space slopea can also be calculated from the slope w of a world-space elevation angle θ
world (i.e., w
=
tan
(
θ
world
)) by using the conversion formula a
=
(
sign
⁡
(
Y
→
y
)
w
s
(
X
→
y
2
+
Y
→
y
2
)
-
w
2
-
s
X
→
y
Y
→
y
)
/
(
s
Y
→
y
2
-
w
2
), where s
=
1
+
w
2 and sign
⁡
(
x
)
=
[
x
≥
0
]
-
[
x
<
0
]. This a is only valid if an equivalent elevation with a positive x component in principal space exists, which is the case if s
(
X
→
y
2
+
Y
→
y
2
)
-
w
2
≥
0 and [
1
,
w
]
·
[
1
-
X
→
y
2
,
X
→
y
]
>
0. This world-space elevation conversion may be particularly useful when used together with the next two principal-space planners.


## 6.3. Hitting the Target Given the Initial Slope

This subsection is about finding the timet
r at which a projectile will hit position [
r
x
,
r
y
] while being launched at slope a. Planning this way may be useful when there is control over the projectile’s initial speed but not over its direction (e.g., for some weapon mounted on a fixed rig).This problem is actually simply a special case of the previous planner, whereb
=
0, meaning that the problem is equivalent to finding the trajectory that touches the line y
=
a
x. After substituting b in (24) and applying some basic algebra, the solution may be written more compactly as follows:
(25)
t
r
=
s
+
s
(
s
+
2
v
∞
)
k
v
∞
,
where s
=
(
1
/
2
)
k
(
a
r
x
-
r
y
). Obviously, the slope a has to be steeper than the line from the origin to the target position (i.e., a
>
r
y
/
r
x) for a solution to exist, which is the case for all but one target position in Figure 12.Figure 12
Planning trajectories through six different targets, all starting with the slope 3/4.


## 6.4. Hitting the Target Given the Target Slope

Similarly, it is possible to hit[
r
x
,
r
y
] given the exact slope a at the target position. This type of planning allows for exact control over the angle at which a target is hit. Again, this can be seen as a special case of the planner from Section 6.2, using r
y
-
a
r
x as the value for b. When substituted and simplified, this results in the following more direct formula:
(26)
t
r
=
r
y
-
a
r
x
k
v
∞
.For this particular planner, the slopea needs to be less than the slope of the line from the origin to the target position (i.e., a
<
r
y
/
r
x). Consequently, only five of the six targets have a valid trajectory in Figure 13.Figure 13
Planning trajectories through six different targets, all ending at the slope −3/4.


## 6.5. Hitting the Target Given the Arc Height or “Curviness”

The timet
r can also be calculated for a target position [
r
x
,
r
y
] and an “arc height” b. Here, b is defined as the maximum difference in the y direction between the trajectory and the straight line from origin to target. Equivalently, b may be interpreted as the height in principal space of the smallest parallelogram containing the whole trajectory, as shown in Figure 14. This problem is another special case of the “line touching” problem and can be solved by substituting r
y
/
r
x for a in (24). After applying some basic algebra, the resulting formula may be written as follows:
(27)
t
r
=
b
v
∞
+
2
b
k
v
∞
.Figure 14
Planning trajectories through six differentRs, all with an arc height of 50.This function is particularly intuitive to plan with when(28)
b
=
h
r
x
2
+
r
y
2
,
where h defines the “curviness” of the trajectory. For example, h
=
0.01 always leads to a low arc for any target position, while h
=
0.5 always leads to a fairly high arc.


## 6.6. Hitting the Target with (Almost) Minimal Effort

A trajectory can also be planned to hit a target at[
r
x
,
r
y
] with the smallest initial speed (and thus the least amount of energy) possible. Planning a trajectory this way requires finding the positive time-to-target t
r which solves ∂
(
v
0
,
x
2
+
v
0
,
y
2
)
/
∂
t
r
=
0, where v
0
,
x and v
0
,
y are defined by (21) and (22), respectively. This equation can be expanded into the following form:
(29)
(
k
2
v
∞
2
)
t
r
4
+
(
k
2
v
∞
r
y
)
t
r
3
-
k
(
r
x
2
+
r
y
2
)
t
r
-
(
r
x
2
+
r
y
2
)
=
0
.Like all quartic equations, solving this in a closed form is possible but difficult to do robustly [10]. In practice, quartic equations are typically solved for any or all of their roots by generic iterative root solvers [11]. But by exploiting domain-specific knowledge, it is also possible to implement a specialized iterative solver for (29) that is guaranteed to efficiently converge to the right root directly. One possible implementation is presented as a C++ function calledGetTimeToTargetRWithMinimalInitialSpeed() in Algorithm 1. There, the equation is first transformed into an equivalent but more well-behaved strictly convex quartic function for which a conservative initial guess for t (or rather, u) is calculated, which is then refined using multiple (but typically less than a total of six) modified or normal conservative Newton’s method iterations. This implementation has been carefully crafted with both robustness and efficiency in mind. As with any implementation, numeric precision can become an issue when using extreme values, but results for practical ranges are typically within a fewfloat epsilons of the exact value. See the comments in the implementation itself for more details.Alternatively, when only a rough approximation of the minimal effort solution is needed, the simpler “curviness” planner from (27) and (28) could be used with h
=
1
/
4. This approximation is fairly accurate for larger values of v
∞, as it actually converges perfectly to the exact solution when v
∞ goes to infinity. But for high friction scenarios, the difference between the exact method and the approximation becomes quite noticeable. To get an idea of the size of the error for a medium friction scenario, the trajectories resulting from the exact method and the approximation are shown side by side in Figure 15 for the case of g
=
9.81m/s2 and v
∞
=
25m/s.Figure 15
Planning trajectories through six differentRs while minimizing the initial velocity using the exact method (in black) and the approximate h
=
1
/
4 method (in red) in a medium-fiction scenario.


## 6.7. Hitting a Target Given the Initial Speed

The last planning algorithm covered in this paper solves for the timet
r required to hit the target [
r
x
,
r
y
] given the projectile’s exact initial speed s. The solution is found by solving v
0
,
x
2
+
v
0
,
y
2
=
s
2 for the positive time t
r to position [
r
x
,
r
y
], where v
0
,
x and v
0
,
y are once again defined by (21) and (22), respectively. This particular equation may be expanded into the following quartic function:
(30)
(
k
2
v
∞
2
)
t
r
4
+
(
2
k
2
v
∞
r
y
)
t
r
3
+
(
k
2
r
x
2
+
k
2
r
y
2
+
2
k
v
∞
r
y
-
s
2
)
t
r
2
+
2
k
(
r
x
2
+
r
y
2
)
t
r
+
(
r
x
2
+
r
y
2
)
=
0
.Note that whens is smaller than the minimal initial speed to hit [
r
x
,
r
y
] (i.e., the sought solution of (29)), the problem has no valid solutions. But when s is larger than that, it will always have exactly two valid solutions. In that case, the smaller of the two t
r values represents the time to hit the target with a low arc, while the larger t
r root is the solution for a high arc. See Figure 16 for an example.Figure 16
Planning trajectories through six differentRs, all starting at 100 m/s. Each of the five reachable targets has a trajectory with a (red) low arc and a (black) high arc.Using a similar approach as in the previous section, a specialized iterative solver can be implemented to solve this particular quartic function. The C++ functionGetTimeToTargetRGivenInitialSpeedS() in Algorithm 2 solves the equation for either the high or low arc root and returns the resulting t
r or returns 0 if no solution exists. This implementation has also been written with robustness, efficiency, and accuracy for a wide range of parameters in mind. Tests showed that the procedure typically requires about six (modified) Newton’s method iterations in total to converge to almost fullfloat precision.


## 7. Future work

As already hinted at in Section4, planning to hit moving targets can sometimes be done using a feedback loop between a target position prediction formula and a planner for a static target, together converging to a solution over multiple iterations. More research is necessary to explore the exact boundary conditions for this convergence to occur or alternatively to look for ways to solve these problems analytically.Additionally, it is likely that the model presented in this paper can also be used for efficient exact collision detection between a trajectory and an arbitrary polygonal mesh by testing the trajectory in principal space against the intersection of the mesh and the principal plane. This algorithm might even be combined with the planners from Sections6.1 and 6.2 to allow for efficient planning of the most optimal trajectory above or below a given polygonal mesh, respectively. These possibilities have not been investigated in depth for this paper but might be covered in future work.


## 8. Conclusion

A novel analytic approximation of ballistic trajectories with air resistance has been presented that was designed to balance physical accuracy and performance in a way that makes sense in the field of computer games.The approximation’s linearity in velocity has been used to define a special principal frame of reference, which makes it possible to always work with these trajectories in a simplified 2D space, even though the original problem can be in 3D with wind coming from any direction. The combined result is that the proposed model is able to produce trajectories that are complex enough to be physically plausible, while keeping the math simple enough to also allow for many different ways of efficient trajectory planning that otherwise might be too impractical for use in computer games.



---

*Source: 463489-2014-01-30.xml*