Add your own alternative version
Stats
217.1K views 3.9K downloads 156 bookmarked
Posted
22 Jun 2007

Comments and Discussions



Hey,
First of all, excellent article and exactly what I was hoping to find when I started researching how to do this for a project I'm working on, so thankyou!
However, your comment about the equation d = (1/2)At^2 (for determining distance after t secs of acceleration at A) was hard to explain and having a mysterious 1/2 in there worried me a bit, so I tried testing it out and I think you might have made an error.
Oh, Secondly, this is my first (ever) visit to codeproject, I don't know its markup, and I couldn't immediately see how to superscript/subscript things, so please bear with the horrible ascii version.
Supposing an acceleration of 3 meters per second per second:
At time t0, velocity is 0, distance is 0
At time t1, velocity is 3, distance is t0+3=3
At time t2, velocity is 6, distance is t1+6=9
At time t3, velocity is 9, distance is t2+9=18
At time t4, velocity is 12, distance is t3+12=30
Trying your Guidance System UI, I configured the following:
Target acceleration (0,0,0)
Target velocity (0,0,0)
Target position (30,0,0)
Source velocity (0,0,0)
Source position (0,0,0)
Projectile acceleration 3
Projectile velocity 0
Projectile pos 0 (not sure what this is?)
Projectile burn time 999
This gives me a stage 1 intersection at 4.47214, surely the correct answer is 5?
Try this instead:
at time t, we know that d = A * (1+2+..+t)
to calculate (1+2+..+t), break the set of additives into pairs:
highest number (t) + lowest number (1) = t+1
next highest (t1) + next lowest (2) = t+1
... and so on
if t is an even number, then we have t/2 pairs of (t+1), so sum is (t+1)*t/2
if t is an odd number, then we have (t1)/2 pairs of (t+1), and a "middle number"
middle number = (t+1)/2
Sum of pairs = (t+1)*(t1)/2
Sum of pairs = (t^2  t + t  1)/2 (expanding)
Sum of pairs = (t^2  1)/2 (simplifying)
now adding sum of pairs and middle number together:
total = (t^2  1)/2 + (t+1)/2
total = (t^2 + t)/2
total = t*(t+1)/2
So for both even and odd values of t, we have:
d = A*t*(t+1)/2 (my equation)
d = A*(t+t^2)/2 (my equation, expanded, reordering for clarity)
d = A*(0+t^2)/2 (your equation, reordering for clarity)
thoughts?





I am investigating whether or not your algorithms with work for our application. Our munition simulator uses 6DOF models (done in Simulink outputted as C++ code) created by a vendor but right now they do not guide so well. So I have asked for more generic (kinematic) models that could hit the target (unrealistically if necessary) every time.
Your algorithms look like the ticket. We need models for high flying guided bomb units (aka, smart bombs) and Hellfire missiles. The main difference between both is the hellfires have initial thrust. After the thrust period it falls with gravity (and drag but we assume its 0 since the models are simpler than 6DOF).
Can your algorithm be modified for the gravity only GBUs? There is initial velocity for both. Targets can be fixed or moving.
Steve
modified 9Oct13 11:20am.





Should it not be
( Pam * Rbt + Pvm ) * t + ( (1/2) * Pam * Rbt^2 + Pvm * Rbt + Ppm )
?





It has been so long since I first wrote this that I have stopped frequently checking for replies. My apologies there.
I presume you are referencing the lines
double ssC = tC  pam * pam * rbt * rbt  2.0 * pam * pvm * rbt  pvm * pvm;
double ssD = tD + pam * pam * rbt * rbt * rbt + pam * pvm * rbt * rbt  2.0 * (pam * rbt + pvm) * ppm;
double ssE = tE  0.25 * pam * pam * rbt * rbt * rbt * rbt + pam * ppm * rbt * rbt  ppm * ppm;
in conjunction with
List<double> ssroots = rootsOf(tA, tB, ssC, ssD, ssE);
as roots correspond to times.
This project was undertaken some years ago, but let's see if I can't still throw down some numbers.
Abstractly, we wish to compute T(t) = P(t), that is, the target and projectile are at the same point in space (the functions' output, T and P) at the same time (the functions' input, t), the point being that if our target and projectile are at the same location at the same time, we've "won."
So, on with the showP(t) for postrbt times is derived as follows:
In plain English, we are the distance that we go while burning all of our fuel, plus our output speed after that burn times the remaining time afterwords, away from where we started.
Or:
P =
(0.5*pam*rbt*rbt + pvm*rbt + ppm)
+
((pam*rbt + pvm) * (trbt))
the first piece is the entire distance covered by the first stage, and the second piece is how far one drifts during the second stage. Technically speaking, over in math world, it should be noted that this function's domain is limited for times of t >= rbt, but we already kinda said that when we described what it is used forthe second stage, and thus why this equation is not shown on the left side of that vertical gray line in the time graph gui.
expanding everything out we get P = 0.5*pam*rbt*rbt + pvm*rbt + ppm + pam*rbt*t  pam*rbt*rbt + pvm*t  pvm*rbt and that can be condensed to P = 0.5*pam*rbt*rbt + ppm + pam*rbt*t + pvm*t by combining like terms
Since the T(t) equation is the square root of a sum of squares (a la Pythagoras), we square both sides of our original T(t) = P(t) to eliminate the radical that appears in T(t), and that leaves us with
P(t)^2 = (0.5*pam*rbt*rbt + ppm + pam*rbt*t + pvm*t) * (0.5*pam*rbt*rbt + ppm + pam*rbt*t + pvm*t)
Now it's just a matter of elbow grease. The proof is *not* trivial, so it won't be left as an exercise to the reader. Multiplying each term by each other term, and using the shortcut that any term multiplied by itself will be itself squared, and any A term multiplied by an B term will appear as both A*B and B*A so for the sake of simplicity 2*A*B can be used insteadI've done a bit of coefficient math in my head, but nothing too fancycertainly nothing beyond 1 = 2*(0.5) or (+)0.25 = (0.5)*(0.5)
P(t)^2 =
0.25*pam*pam*rbt*rbt*rbt*rbt  pam*rbt*rbt*ppm  pam*rbt*rbt*pam*rbt*t  pam*rbt*rbt*pvm*t
+
ppm*ppm + 2*ppm*pam*rbt*t + 2*ppm*pvm*t
+
pam*rbt*t*pam*rbt*t + 2*pam*rbt*t*pvm*t
+
pvm*t*pvm*t
collecting all degrees of t together...
P(t)^2 =
(pam*pam*rbt*rbt + 2*pam*pvm + pvm*pvm)*t*t
+
(2*pam*ppm*rbt + 2*pvm*ppm  pam*pam*rbt*rbt*rbt  pam*pvm*rbt*rbt)*t
+
(0.25*pam*pam*rbt*rbt*rbt*rbt  pam*ppm*rbt*rbt + ppm*ppm)
It's a miracleI check my work years later and I don't make a math error on my first try. When you flip all of the signs (since this function is being subtracted from T(t)^2), you will clearly see that this corresponds to the projectilerelevant parts of ssC, ssD, and ssE, as ultimately we are solving for 0 = T(t)^2 + P(t)^2, in the most verbose terms.
And my math is vindicated in the little minigame I provided; missiles shot without enough fuel to impact in stage one, but which will gain enough speed to hit in stage two, do in fact hit.
If you'd like to share how you arrived at your equation perhaps I can shed some light on what is going on. And again sorry for not getting to this soonerI am, after all, so busy with college these days. It is still possible that there may be an error in the article itself, even if the code is not flawed, so if you're picking up on something there, be specific about it and I'll see if it corresponds to my code as it should. If any error exists in the text I will have it corrected as soon as possible.
Gatsby
modified on Wednesday, September 1, 2010 6:53 AM





http://www.flickr.com/photos/27513972@N06/2564214779/





I ported both your code and the code from the UnReal Wiki (http://wiki.beyondunreal.com/Legacy:Projectile_Aiming ) to c++ and ran a comparison. The test is exactly this:
Shooter is a CIWIS on a Iowa class battleship moving at about 10 knots.
Target is a Silkworm missile moving at roughly 500 knots.
At ranges less than 1000 yards both solutions are almost identical. However as the range increases the Unreal Wiki solution consistently gets closer to the target. At > 10000 yards I can't even see your solutions projectiles they're so far away. I can do a webex with you if you want or send you screen shots.
I'm not a mathematician so I don't understand the difference but is your solution algabraic and the other another type(geometric?). The other is using a lot of Dot and Cross products is the only difference I can see.
Any way your solution is much more stable, the UnrealWiki solution produces a lot of erroneous trajectories (but the valid one always hit smack on target). It's very strange and I'd really like to figure it out because your solution is much more stable and generates fewer errors.
John





I most certainly do wish to troubleshoot this with you, but I feel I need a little more information. Any sort of code you can share with me? I'm particularly interested in what you're passing in to my functions, and what sort of environmental factors may be at work (is there air resistance, is there any sort of gravity, etc)
Please reply back and I'll be very happy to help.
EDIT: ah yes, Also, I ported my code to C++ myself recently, allow me to provide that:
MGS.h
#include <cmath>
#include <vector>
struct Point3D
{
double x;
double y;
double z;
Point3D();
Point3D(double X, double Y, double Z);
};
struct Projectile
{
double a;
double v;
double p;
Projectile();
Projectile(double A, double V, double P);
};
struct Shooter
{
Point3D v;
Point3D p;
Shooter();
Shooter(Point3D V, Point3D P);
};
struct Target
{
Point3D a;
Point3D v;
Point3D p;
Target();
Target(Point3D A, Point3D V, Point3D P);
};
Point3D CalculateAcceleration(Point3D PAtTMinus2, Point3D PAtTMinus1, Point3D PAtTMinus0);
Point3D CalculateVelocity(Point3D PAtTMinus2, Point3D PAtTMinus1, Point3D PAtTMinus0);
bool Intercept(Target T, Shooter S, Projectile P, int rbt, double tolerance, std::vector<Point3D>& vectors);
MGS.cpp
#include <algorithm>
#include <limits>
#include "MGS.h"
#define EPSILON 0.00000000001
using namespace std;
static inline bool isnan(double value)
{
return value != value;
}
static inline bool isinf(double value)
{
return value == std::numeric_limits<double>::infinity();
}
Point3D::Point3D()
{
x = 0;
y = 0;
z = 0;
}
Point3D::Point3D(double X, double Y, double Z)
{
x = X;
y = Y;
z = Z;
}
Projectile::Projectile()
{
a = 0;
v = 0;
p = 0;
}
Projectile::Projectile(double A, double V, double P)
{
a = A;
v = V;
p = P;
}
Shooter::Shooter()
{
v = Point3D();
p = Point3D();
}
Shooter::Shooter(Point3D V, Point3D P)
{
v = V;
p = P;
}
Target::Target()
{
a = Point3D();
v = Point3D();
p = Point3D();
}
Target::Target(Point3D A, Point3D V, Point3D P)
{
a = A;
v = V;
p = P;
}
Point3D CalculateAcceleration(Point3D PAtTMinus2, Point3D PAtTMinus1, Point3D PAtTMinus0)
{
return Point3D(PAtTMinus2.x  2.0 * PAtTMinus1.x + PAtTMinus0.x,
PAtTMinus2.y  2.0 * PAtTMinus1.y + PAtTMinus0.y,
PAtTMinus2.z  2.0 * PAtTMinus1.z + PAtTMinus0.z);
}
Point3D CalculateVelocity(Point3D PAtTMinus2, Point3D PAtTMinus1, Point3D PAtTMinus0)
{
return Point3D(PAtTMinus2.x / 2.0  2.0 * PAtTMinus1.x + 3.0 * PAtTMinus0.x / 2.0,
PAtTMinus2.y / 2.0  2.0 * PAtTMinus1.y + 3.0 * PAtTMinus0.y / 2.0,
PAtTMinus2.z / 2.0  2.0 * PAtTMinus1.z + 3.0 * PAtTMinus0.z / 2.0);
}
static vector<double> rootsOf(double A, double B, double C, double D, double E)
{
vector<double> roots = vector<double>();
if (A == 0.0)
{
if (B == 0.0)
{
if (C == 0.0)
{
if (D == 0.0)
{
if (E == 0.0)
roots.push_back(0.0);
}
else
roots.push_back(D / E);
}
else
{
roots.push_back((D + sqrt(D * D  4 * C * E)) / (2.0 * C));
roots.push_back((D  sqrt(D * D  4 * C * E)) / (2.0 * C));
}
}
else
{
C /= B;
D /= B;
E /= B;
double F = (3.0 * D  C * C) / 3.0;
double G = (2.0 * C * C * C  9.0 * C * D + 27.0 * E) / 27.0;
double H = (G * G) / 4.0 + (F * F * F) / 27.0;
if (H > 0)
{
double intermediate = G / 2.0 + sqrt(H);
double m = intermediate < 0.0 ? pow(intermediate, 1.0 / 3.0) : pow(intermediate, 1.0 / 3.0);
intermediate = 2.0 * sqrt(H);
double n = intermediate < 0.0 ? pow(intermediate, 1.0 / 3.0) : pow(intermediate, 1.0 / 3.0);
roots.push_back(m + n  C / 3.0);
}
else
{
double intermediate = sqrt(G * G / 4.0  H);
double rc = intermediate < 0.0 ? pow(intermediate, 1.0 / 3.0) : pow(intermediate, 1.0 / 3.0);
double theta = acos(G / (2.0 * intermediate)) / 3.0;
roots.push_back(2.0 * rc * cos(theta)  C / 3.0);
roots.push_back(rc * (cos(theta) + sqrt(3.0) * sin(theta))  C / 3.0);
roots.push_back(rc * (cos(theta)  sqrt(3.0) * sin(theta))  C / 3.0);
}
if (F + G + H == 0.0)
{
double intermediate = E < 0.0 ? pow(E, 1.0 / 3.0) : pow(E, 1.0 / 3.0);
roots.clear();
roots.push_back(intermediate);
roots.push_back(intermediate);
roots.push_back(intermediate);
}
}
}
else
{
B /= A;
C /= A;
D /= A;
E /= A;
double F = C  (3.0 * B * B) / 8.0;
double G = D + B * B * B / 8.0  (B * C) / 2.0;
double H = E  3.0 * B * B * B * B / 256.0 + B * B * C / 16.0  B * D / 4.0;
double b = F / 2.0;
double c = (F * F  4.0 * H) / 16.0;
double d = (G * G) / 64.0;
double f = (3.0 * c  b * b) / 3.0;
double g = (2.0 * b * b * b  9.0 * b * c + 27.0 * d) / 27.0;
double h = (g * g) / 4.0 + (f * f * f) / 27.0;
double y1;
double y2r;
double y2i;
double y3r;
double y3i;
if (h > 0.0)
{
double intermediate = g / 2.0 + sqrt(h);
double m = intermediate < 0.0 ? pow(intermediate, 1.0 / 3.0) : pow(intermediate, 1.0 / 3.0);
intermediate = 2.0 * sqrt(h);
double n = intermediate < 0.0 ? pow(intermediate, 1.0 / 3.0) : pow(intermediate, 1.0 / 3.0);
y1 = m + n  b / 3.0;
y2r = (m + n) / 2.0  b / 3.0;
y2i = ((m  n) / 2.0) * sqrt(3.0);
y3r = (m + n) / 2.0  b / 3.0;
y3i = ((m  n) / 2.0) * sqrt(3.0);
}
else
{
double intermediate = sqrt((g * g / 4.0  h));
double rc = intermediate < 0.0 ? pow(intermediate, 1.0 / 3.0) : pow(intermediate, 1.0 / 3.0);
double theta = acos((g / (2.0 * intermediate))) / 3.0;
y1 = 2.0 * rc * cos(theta)  b / 3.0;
y2r = rc * (cos(theta) + sqrt(3.0) * sin(theta))  b / 3.0;
y2i = 0.0;
y3r = rc * (cos(theta)  sqrt(3.0) * sin(theta))  b / 3.0;
y3i = 0.0;
}
if (f + g + h == 0.0)
{
double intermediate = d < 0.0 ? pow(d, 1.0 / 3.0) : pow(d, 1.0 / 3.0);
y1 = intermediate;
y2r = intermediate;
y2i = 0.0;
y3r = intermediate;
y3i = 0.0;
}
double p;
double q;
if (h <= 0.0)
{
int zeroCheck = 0;
vector<double> cubicRoots(3);
cubicRoots[0] = y1;
cubicRoots[1] = y2r;
cubicRoots[2] = y3r;
sort(cubicRoots.begin(), cubicRoots.begin() + 2);
p = sqrt(cubicRoots[1]);
q = sqrt(cubicRoots[2]);
if (abs(y1) < EPSILON)
{
p = sqrt(y2r);
q = sqrt(y3r);
zeroCheck = 1;
}
if (abs(y2r) < EPSILON)
{
p = sqrt(y1);
q = sqrt(y3r);
zeroCheck += 2;
}
if (abs(y3r) < EPSILON)
{
p = sqrt(y1);
q = sqrt(y2r);
zeroCheck += 4;
}
switch (zeroCheck)
{
case (3):
p = sqrt(y3r);
break;
case (5):
p = sqrt(y2r);
break;
case (6):
p = sqrt(y1);
break;
}
if (y1  EPSILON < 0.0  y2r  EPSILON < 0.0  y3r  EPSILON < 0.0)
{
if (E == 0.0)
roots.push_back(0.0);
}
else
{
double r;
if (zeroCheck < 5)
r = G / (8.0 * p * q);
else
r = 0.0;
double s = B / 4.0;
roots.push_back(p + q + r  s);
roots.push_back(p  q  r  s);
roots.push_back(p + q  r  s);
roots.push_back(p  q + r  s);
}
}
else
{
double r2mod = sqrt(y2r * y2r + y2i * y2i);
double y2mod = sqrt((r2mod  y2r) / 2.0);
double x2mod = y2i / (2.0 * y2mod);
p = x2mod + y2mod;
double r3mod = sqrt(y3r * y3r + y3i * y3i);
double y3mod = sqrt((r3mod  y3r) / 2.0);
double x3mod = y3i / (2.0 * y3mod);
q = x3mod + y3mod;
double r = G / (8.0 * (x2mod * x3mod + y2mod * y3mod));
double s = B / 4.0;
roots.push_back(x2mod + x3mod + r  s);
roots.push_back(x2mod  x3mod + r  s);
}
}
for (int i = 0; i != roots.size(); ++i)
if (isinf(roots[i])  isnan(roots[i]))
roots.erase(roots.begin() + i);
sort(roots.begin(), roots.end());
return roots;
}
bool does_not_contain(vector<double> vector, double value)
{
for(int i = 0; i != vector.size(); ++i)
if(vector[i] == value)
return false;
return true;
}
void keysort(vector<double>& keys, vector<double>& items)
{
for(int i = 0; i != keys.size(); ++i)
for(int j = keys.size()  1; j != i; j)
if(keys[j  1] > keys[j])
{
double temp = keys[j];
keys[j] = keys[j  1];
keys[j  1] = temp;
temp = items[j];
items[j] = items[j  1];
items[j  1] = temp;
}
}
bool Intercept(Target T, Shooter S, Projectile P, int rbt, double tolerance, vector<Point3D>& vectors)
{
T.v.x = S.v.x;
T.v.y = S.v.y;
T.v.z = S.v.z;
T.p.x = S.p.x;
T.p.y = S.p.y;
T.p.z = S.p.z;
if (tolerance < 1.0)
tolerance = 1.0;
if (rbt < 0)
rbt = 0;
if (rbt == 0)
P.a = 0.0;
double tA = 0.25 * (T.a.x * T.a.x + T.a.y * T.a.y + T.a.z * T.a.z);
double tB = T.a.x * T.v.x + T.a.y * T.v.y + T.a.z * T.v.z;
double tC = T.v.x * T.v.x + T.a.x * T.p.x + T.v.y * T.v.y + T.a.y * T.p.y + T.v.z * T.v.z + T.a.z * T.p.z;
double tD = 2.0 * (T.v.x * T.p.x + T.v.y * T.p.y + T.v.z * T.p.z);
double tE = T.p.x * T.p.x + T.p.y * T.p.y + T.p.z * T.p.z;
double fsA = tA  0.25 * P.a * P.a;
double fsB = tB  P.a * P.v;
double fsC = tC  P.a * P.p  P.v * P.v;
double fsD = tD  2.0 * P.v * P.p;
double fsE = tE  P.p * P.p;
double ssC = tC  P.a * P.a * rbt * rbt  2.0 * P.a * P.v * rbt  P.v * P.v;
double ssD = tD + P.a * P.a * rbt * rbt * rbt + P.a * P.v * rbt * rbt  2.0 * (P.a * rbt + P.v) * P.p;
double ssE = tE  0.25 * P.a * P.a * rbt * rbt * rbt * rbt + P.a * P.p * rbt * rbt  P.p * P.p;
vector<double> roots = vector<double>();
vector<double> fsroots = rootsOf(fsA, fsB, fsC, fsD, fsE);
vector<double> fsrootsCA = rootsOf(0.0, 4.0 * fsA, 3.0 * fsB, 2.0 * fsC, fsD);
vector<double> ssroots = rootsOf(tA, tB, ssC, ssD, ssE);
vector<double> ssrootsCA = rootsOf(0.0, 4.0 * tA, 3.0 * tB, 2.0 * ssC, ssD);
for(int i = 0; i != fsroots.size(); ++i)
{
double root = fsroots[i];
if(root >= 0.0 && root <= rbt && does_not_contain(roots, root))
roots.push_back(root);
}
for(int i = 0; i != ssroots.size(); ++i)
{
double root = ssroots[i];
if(root > rbt && does_not_contain(roots, root))
roots.push_back(root);
}
for(int i = 0; i != fsrootsCA.size(); ++i)
{
double root = fsrootsCA[i];
if (tolerance >= abs(fsA * root * root * root * root + fsB * root * root * root + fsC * root * root + fsD * root + fsE) && root >= 0.0 && root <= rbt && does_not_contain(roots, root))
roots.push_back(root);
}
for(int i = 0; i != ssrootsCA.size(); ++i)
{
double root = ssrootsCA[i];
if (tolerance >= abs(tA * root * root * root * root + tB * root * root * root + ssC * root * root + ssD * root + ssE) && root > rbt && does_not_contain(roots, root))
roots.push_back(root);
}
sort(roots.begin(), roots.end());
bool able = roots.size() != 0;
if(!able)
{
vector<double> fsattempts = vector<double>();
vector<double> ssattempts = vector<double>();
for(int i = 0; i != fsrootsCA.size(); ++i)
{
double root = fsrootsCA[i];
if (root >= 0.0 && root <= rbt && does_not_contain(fsattempts, root))
fsattempts.push_back(root);
}
fsattempts.push_back(0);
for(int i = 0; i != ssrootsCA.size(); ++i)
{
double root = ssrootsCA[i];
if(root > rbt && does_not_contain(ssattempts, root))
ssattempts.push_back(root);
}
vector<double> keys = vector<double>(fsattempts.size() + ssattempts.size());
vector<double> items = vector<double>(fsattempts.size() + ssattempts.size());
for (int i = 0; i != fsattempts.size(); ++i)
{
keys[i] = abs(fsA * fsattempts[i] * fsattempts[i] * fsattempts[i] * fsattempts[i] + fsB * fsattempts[i] * fsattempts[i] * fsattempts[i] + fsC * fsattempts[i] * fsattempts[i] + tD * fsattempts[i] + tE);
items[i] = fsattempts[i];
}
for (int i = 0; i != ssattempts.size(); ++i)
{
keys[fsattempts.size() + i] = abs(tA * ssattempts[i] * ssattempts[i] * ssattempts[i] * ssattempts[i] + tB * ssattempts[i] * ssattempts[i] * ssattempts[i] + ssC * ssattempts[i] * ssattempts[i] + ssD * ssattempts[i] + ssE);
items[fsattempts.size() + i] = ssattempts[i];
}
keysort(keys, items);
for(int i = 0; i != items.size(); ++i)
roots.push_back(items[i]);
}
vectors = vector<Point3D>(roots.size());
for(int i = 0; i != roots.size(); ++i)
{
Point3D v;
if(P.a == 0.0 && P.v == 0.0)
{
v.x = S.v.x;
v.y = S.v.y;
v.z = S.v.z;
}
else
{
v.x = 0.5 * T.a.x * roots[i] * roots[i] + T.v.x * roots[i] + T.p.x;
v.y = 0.5 * T.a.y * roots[i] * roots[i] + T.v.y * roots[i] + T.p.y;
v.z = 0.5 * T.a.z * roots[i] * roots[i] + T.v.z * roots[i] + T.p.z;
}
double m = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
if(m < EPSILON)
m = 1.0;
vectors[i].x = v.x / m;
vectors[i].y = v.y / m;
vectors[i].z = v.z / m;
}
return able;
}
Testable then by including MGS.h and making some shooter, target, and projectile structs and passing them in.





I'm not using air resistence, I am using gravity. I'm passing that in via the Target's Acceleration which is only downward towards the earth.
Here's my actual code:
std::vector<dvector> Targeting::GetShootVect(
Vector3 StartLoc,
Vector3 EndLoc,
float ProjSpeed,
Vector3 Gravity,
Vector3 launcherVelocity,
Vector3 targetVelocity)
{
double ROT = 400;
std::vector<dvector> vectors;
Intercept(
Gravity.x,Gravity.y,Gravity.z,
targetVelocity.x, targetVelocity.y, targetVelocity.z,
EndLoc.x, EndLoc.y, EndLoc.z,
launcherVelocity.x, launcherVelocity.y, launcherVelocity.z,
StartLoc.x, StartLoc.y, StartLoc.z,
0, ProjSpeed, 0,
500, ROT, vectors);
return vectors;
}





One other thing that you should know about gravity.
The Target is flying at a constant height (it isn't falling). Every time step I'm fixing that so it technically has no acceleration. You would think gravity should be 0 in this case but I find it I set gravity to zero the projectiles pass under the Target at the same height that the Target would be if it were falling. I saw this and then passed in gravity for the Target acceleration and everything worked (except for sometimes the projectiles miss the target which is the reason I think there's some sort of error).
I sort of lumped everything into the same bucket (target and projectile acceleration wise). This may be wrong but I couldn't find anywhere to define the acceleration of the projectile using your code.





One thing I can answer right away is
"I couldn't find anywhere to define the acceleration of the projectile using your code"
 the var "pam" is the one you're looking for. This is the projectile's acceleration vector's magnitude, just as "pvm" is it's initial or "muzzle" velocity, and "ppm" the position offset (were this a person firing a gun, you might think of this as the length of their arm). Just be sure that when you use the pam var that you provide as much Remaining Burn Time as you want  passing in infinity should work if you really don't want the thing to be running out of fuel.
It's difficult for me to find the error give the code I've got, if you'd be so obliged you can email me directly with more of the context in which the code appears and I will be more able to help you [email removed], because as it stands right now I can't figure out where the problem is. I've tested things similar to what you're describing in my own demo apps and everything seems to be working perfectly, so I suspect that it's how the code is being implemented that's where the error is coming in to play (largely my fault, as I didn't give as much code documentation as I probably should have when I released this article  my apologies). If I'm able to see how you're updating the positions of things and where exactly in the updating that the intercept is being calculated, etc etc etc.
Please do continue to reply, as I'm glad you've found my code helpful, and I really would like to get this to work for you. As stated above, I've checked and rechecked my stuff using scenarios such as the ones you've described and everything works fine for me, so it really should just be an implementation bug, which I'd be happy to help track down with you.
[issue resolved]
Gatsby
modified on Friday, June 13, 2008 4:30 PM





I was trying to solve a similar problem but couldn't quite handle the math so I downloaded an example from the UnReal Wiki which did all the work for me. I like your solution better though.
http://wiki.beyondunreal.com/Legacy:Projectile_Aiming
I'm curious what you think of their solution.





First, could you comment your code? Thanks.
Second, I've been trying off and on to modify the intercept code for a matching intercept  ie., one in which, at the end of the intercept, the projectile is going in the same direction, with the same speed, as the target. (Once that is matched, acceleration can be matched outside of the formula and would probably be better matched outside of the formula anyway).
It would be a little too much to ask you to solve it for me, but any help you would want to give would be appreciated.
Thanks.





Hello Jay Gatsby,
First of all, i will say that you wrote a interesting article.
Unfortunately, i do not understand the vector in this line:
Targeting.Intercept(Tracking.aX, Tracking.aY, 0.0, Tracking.vX, Tracking.vY, 0.0, Tracking.pX, Tracking.pY, 0.0, this.vX, this.vY, 0.0, this.pX, this.pY, 0.0, this.aM, 0.0, 0.0, rbt, ROT, out vectors, out tui);
Is that the aiming(target) point (where the missile will go to)?
If not is there a calculation in your system, that can give me these point in level coordinates?
thx hulk





There is no aiming(target) point, as you have worded it, but there is a static Tracking class which monitors the position of the target every timestep, and, given this history, calculates what acceleration, velocity, and position the target must have to produce the last three locations it was given. I'm careful to avoid calling these computations "estimates," as in my system they are not (due to the fact that objects do in fact behave in a second derivative fashion), but if the target behaved in some unexpected manor, these would be only a best guess (as in times when the target is evading, or in scenarios where the target does not behave in a .5At^2 + Vt + P way).
In a nutshell, We're passing in what appears to be the target's acceleration, velocity, and position, along with all the other things we know. If the target behaves in a constant acceleration fashion, these will be spot on. The reason I did this was to demonstrate that the guidance system does not need the target to be faxing it's movement parameters (like it's current acceleration vector) to the missile in order for the system to function. Rather, I aimed to prove that the missile could figure out what the target is doing if the target isn't so nice as to tell it directly. This is far more realistic (and ergo, not trivial), as my algorithm does not assume that the target related information is directly known. though it is perfectly capable of mathematically solving for what it needs, as it would need to in a practical application.
Gatsby





Hallo Jay Gatsby,
thanks for replying, I work on a small space shooter game similarly yours (but with DX9).My intention was to combine your guidance system with a pathfinding system like A*. But if there is no point wich i can give the pathfinding system this does not work i together.
Hulk





After some thought, you may be able to modify it to give you what you're looking for.
you may recall this block of code from the bottom of the intercept function in Targeting.cs
for (int i = 0; i != roots.Count; i++)
{
double x;
double y;
double z;
if (pam == 0.0 && pvm == 0.0)
{
x = svx;
y = svy;
z = svz;
}
else
{
x = 0.5 * tax * roots[i] * roots[i] + tvx * roots[i] + tpx;
y = 0.5 * tay * roots[i] * roots[i] + tvy * roots[i] + tpy;
z = 0.5 * taz * roots[i] * roots[i] + tvz * roots[i] + tpz;
}
double m = Math.Sqrt((x * x + y * y + z * z));
if (m < double.Epsilon)
m = 1.0;
vectors[i, 0] = x / m;
vectors[i, 1] = y / m;
vectors[i, 2] = z / m;
}
What that's doing is finding out where the point of impact will be, and giving a unit vector in that direction (always do the else code in your case, the if is there merely to point the missile prograde when out of fuel). Note that x, y, and z are *offsets* from the shooter, so you'd have to add those back in if you want global coords. What you'd need to do is pass those offset, ununormalized raw coordinates. I threw them out in the stock intercept method because a unit vector giving direction is far more useful to someone scaling the initial launch vectors higher on the call stack. But yea, the x, y, and z as they come raw from the else block here, added to the spx, spy, and spz values respectively, should be what you're looking for, if I understood your question correctly.
Gatsby







Could you please provide some documentation as to what the arguments your expecting for the functions are? For example, the code you posted in the WW2 message:
Calculate(double Sx, double Sy, double Sh, double Sv, double Tx, double Ty, double Th, double Tv,double Pv, bool InDegrees)
I have no ideal what I'm suppose to pass for Sx, Sy, Sh ...etc.
public static double Calculate(double Sx, double Sy, double Sh, double Sv, double Tx, double Ty, double Th, double Tv, double Pv, bool InDegrees)
What I'm trying to accomplish is given two points with angle and speed, what's the equation to calculate at which point the two points will intersect. I need to make a boolean function that will calculate this and return true if they will intersect, and false if not. Then handle false by adjusting the direction in which point two travels to accommodate an interception with point one.





Good news, I can definitely help you!
I will start by identifying the intent of the parameters you have listed, but bear in mind that they are for firing on a target, not what you have explained what you need to do. Once I've done that I'll detail what you should do in your case.
Sx: the x coordinate of the shooter
Sy: the y coordinate of the shooter
Sh: the heading angle of the shooter
Sv: the speed of the shooter (along the bearing provided by Sh)
Tx, Ty, Th, and Tv are the same as the above four, but for the target, not the shooter
Pv: the speed of the projectile (we don't know it's heading, that's what we'll solve for)
InDegrees: true indicates that Sh and Th (and the returned Ph value) will be in degrees (i.e. 180), while false indicates that these values should be interpreted as radians (i.e. 3.14159...etc.).
What you need to do is pass in obj1.x for Sx, obj1.y for Sy, 0 for Sh, 0 for Sv, obj2.x for Tx, obj2.y for Ty, obj2.h for Th, obj2.v for Tv and obj1.v for Pv (and true or false for InDegrees, depending on weather you're using degrees or not for your angle measures). If an exception is thrown it means that not only will an intercept never occur, there's no way that you can adjust the heading of obj1 ot make an intercept occur. If, however, no exception is thrown, the return value you get will indicate the heading that the first object would need to have in order to intercept, and if the first object's heading is very close to this it will hit second object (i.e. the difference in heading is less than some predefined small constant, 0.000000001 being a good choice. Use AbsoluteValue(h1h2) < smallConstant for this determination, true == headings are really close, false == they're not. Do NOT use zero, rounding errors will plague youthe result will always be false!) You said that you wanted to know the heading that would be needed to intercept if no intercept would occur, well look no further, that Ph return value is what you are looking for. (Please remember to credit me with coming up with this algorithm.)
Let me know if that helps. If you are still having difficulties, please let me know, and provide what you've got so farI'll be glad to help!
Thanks for writing!
Gatsby





Thanks alot, worked like a charm , what I'm working on is a homing missle in 2D for a video game I'm trying to make. Heres what I got for the homing missle so far. Keep in mind that this is sloppy coding and my goal right now was just to get it working. My next step will be to go through and get the rest working then fine tune it. Plus, this only works for the 1st and 2nd quadrant right now. The missle will actually curve towards the target on a fixed radius until the point of intercept is true, so you get more of a realistic effect as the missle travels. This is using the TorqueX game engine and XNA.
public float Intercept(float Sx, float Sy, float Sh, float Sv, float Tx, float Ty,<br />
float Th, float Tv, float Pv, bool InDegrees)<br />
{<br />
if(InDegrees)<br />
{<br />
Sh = Sh * (float)(Math.PI / 180.0);<br />
Th = Th * (float)(Math.PI / 180.0);<br />
}<br />
float theta = (float)Math.Atan2(Tx  Sx, Ty  Sy);<br />
float shOff = Sh  theta;<br />
float thOff = Th  theta;<br />
float svxOff = (float)(Sv * Math.Sin(shOff));<br />
float svyOff = (float)(Sv * Math.Cos(shOff));<br />
float tvxOff = (float)(Tv * Math.Sin(thOff));<br />
float tvyOff = (float)(Tv * Math.Cos(thOff));<br />
float phOff = (float)Math.Asin((tvxOff  svxOff) / Pv);<br />
float crosscheck = (float)Math.Cos(phOff) * Pv;<br />
if(float.IsNaN(phOff)  tvyOff  svyOff >= crosscheck)<br />
return 1;<br />
return InDegrees ? (float)((phOff + theta) * 180.0 / Math.PI) : phOff + theta;<br />
}<br />
<br />
<br />
if((Game.Instance.Missile.Position.X + 2.5f >= Game.Instance.Enemy.Position.X) &&<br />
(Game.Instance.Missile.Position.X  2.5f <= Game.Instance.Enemy.Position.X) &&<br />
(Game.Instance.Missile.Position.Y + 2.5f >= Game.Instance.Enemy.Position.Y) &&<br />
(Game.Instance.Missile.Position.Y  2.5f <= Game.Instance.Enemy.Position.Y))<br />
{<br />
Vector2 vPos = new Vector2(0.0f, 0.0f);<br />
Game.Instance.Missile.SetPosition(vPos, false);<br />
Game.Instance.Missile.SetRotation(180.0f, false);<br />
Globals.setOrigin = true;<br />
}<br />
<br />
float fIntr = Intercept(Game.Instance.Missile.Position.X, Game.Instance.Missile.Position.Y,<br />
Game.Instance.Missile.Rotation, 0.25f, Game.Instance.Enemy.Position.X,<br />
Game.Instance.Enemy.Position.Y, Game.Instance.Enemy.Rotation, 0.1f, 0.25f, true);<br />
<br />
float angle = T2DVectorUtil.AngleFromTarget(Game.Instance.Enemy.Position,<br />
Game.Instance.Missile.Position);<br />
<br />
if((angle > 180.0f) && (angle < 275.0f))
{<br />
if((Globals.fAdjAngle > 275.0f)  (Globals.fAdjAngle < 180.0f))<br />
Globals.fAdjAngle = 275.0f;<br />
if(Globals.setOrigin)<br />
{<br />
Globals.fOrigin = Game.Instance.Missile.Position.X  15.0f;<br />
Globals.setOrigin = false;<br />
}<br />
if(fIntr > angle + 3.5f)<br />
{<br />
Vector2 vPos = new Vector2(Globals.fOrigin * (float)Math.Sin(Globals.DegreesToRadians(Globals.fAdjAngle)),<br />
Globals.fOrigin * (float)Math.Cos(Globals.DegreesToRadians(Globals.fAdjAngle)));<br />
Game.Instance.Missile.SetPosition(vPos, false);<br />
Game.Instance.Missile.SetRotation(angle + 180.0f, false);<br />
}<br />
else<br />
{<br />
Vector2 vPos = new Vector2(Game.Instance.Missile.Position.X + (0.25f * (float)Math.Sin(Globals.DegreesToRadians(Game.Instance.Missile.Rotation))),<br />
Game.Instance.Missile.Position.Y  (0.25f * (float)Math.Cos(Globals.DegreesToRadians(Game.Instance.Missile.Rotation))));<br />
Game.Instance.Missile.SetPosition(vPos, false);<br />
Game.Instance.Missile.SetRotation(angle + 180.0f, false);<br />
}<br />
<br />
Globals.fAdjAngle = 0.25f;<br />
}<br />
else if((angle > 90.0f) && (angle < 180.0f))
{<br />
if((Globals.fAdjAngle > 180.0f)  (Globals.fAdjAngle < 90.0f))<br />
Globals.fAdjAngle = 90.0f;<br />
if(Globals.setOrigin)<br />
{<br />
Globals.fOrigin = Game.Instance.Missile.Position.X + 15.0f;<br />
Globals.setOrigin = false;<br />
}<br />
if(fIntr > 0)<br />
{<br />
Vector2 vPos = new Vector2(Globals.fOrigin * (float)Math.Sin(Globals.DegreesToRadians(Globals.fAdjAngle)),<br />
Globals.fOrigin * (float)Math.Cos(Globals.DegreesToRadians(Globals.fAdjAngle)));<br />
Game.Instance.Missile.SetPosition(vPos, false);<br />
Game.Instance.Missile.SetRotation(angle + 180.0f, false);<br />
}<br />
else<br />
{<br />
Vector2 vPos = new Vector2(Game.Instance.Missile.Position.X + (0.25f * (float)Math.Sin(Globals.DegreesToRadians(Game.Instance.Missile.Rotation))),<br />
Game.Instance.Missile.Position.Y  (0.25f * (float)Math.Cos(Globals.DegreesToRadians(Game.Instance.Missile.Rotation))));<br />
Game.Instance.Missile.SetPosition(vPos, false);<br />
Game.Instance.Missile.SetRotation(angle + 180.0f, false);<br />
}<br />
<br />
Globals.fAdjAngle += 0.25f;<br />
}code><br />
<br />
<br />
 modified at 20:02 Monday 27th August, 2007





YES! *does a victory dance* I helped somebody! Glad to know it worked for you!
<whisper>You wouldn't mind putting "//Intercept algorithm developed by Christopher Kampas, © 2007, all rights reserved" above "public float Intercept(float Sx, float Sy, float Sh, float Sv, float Tx, float Ty, float Th, float Tv, float Pv, bool InDegrees)" would you? *blush* It would be an honor to "see my name in lights," so to speak... you know what I mean?</whisper>
Gatsby





Lol, I can do that, thanks again!





I appreciate it, and thank you!
Over and out.
Gatsby





Jay,
You are missing a point here. Especially, if you ignore impossible cases.
Try solving the missile 'torpedo' problem (in two dimensions) for all possible vessel and weapon speeds using quadratitics or trig. It won't work.
There is a far simpler logical solution, that was sorted out in WWII to produce the mechanical torpedo aiming calculator known as the 'fuit machine'.
No 1GHz processors in those days.
Roland







General News Suggestion Question Bug Answer Joke Praise Rant Admin Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

