Troubleshooting Professional Magazine
|
Simulation Hello World |
| Steve Litt is the author of Troubleshooting Techniques of the Successful Technologist and Rapid Learning: Secret Weapon of the Successful Technologist. |
|
|
It doesn't matter. By reviewing all the articles in this TPM issue, you'll learn the basics of simple trajectory simulation. From there you can improve the design. And from there, given time and desire, you can go as far into simulation as you want. This magazine will get you through the coding of a 3D Vector class, a Body class acted upon by both gravity and wind, and a simple if slightly buggy simulation based aiming program.
The original intent had been to include a missile launcher and missile aimer capable of accuracy half way around the world. Unfortunately, the necessary spherical coordinate math proved so time consuming that it couldn't be completed in time. But perhaps the first issue of the new Linux Productivity Magazine will feature a missile launcher and missile guider simulation. It all depends on time.
Yes, Troubleshooting Professional Magazine will become two separate magazines. Your votes overwhelmingly confirmed my suspicion that all readers would be much better served by two magazines. The details are in the next article.
So kick back, relax, and enjoy this last ever Troubleshooting/Linux combination magazine. No matter what your interest, I think you'll enjoy it.
The vote was unanimous. Everyone wants separate magazines. Most also said that it would be really cool if the troubleshooting mag could be monthly.
So the introductory issue of Linux Productivity Magazine will come out in August 2002. It will be a monthly magazine. Troubleshooting Professional Magazine will return to its roots of "all troubleshooting all the time", and for the time being will be quarterly because I have limited time. Each year there will be a January, April, July and October issue. I'll watch my web stats closely, and if Troubleshooting Professional gets a lot of traffic, I'll find a way to promote it to a monthly mag. If you like the Troubleshooting content in Troubleshooting Professional Magazine, you can vote for it to become monthly just by visiting it.
Linux Productivity Magazine will be formatted and themed very similarly to Troubleshooting Professional Magazine. Obviously there won't be a Linux Log article now that the entire magazine is devoted to Linux and Open Source, but at least for a while there will still be a Life After Windows column.
Loyal readers, thank you for 5 1/2 years of faithful support, and for enthusiastically embracing Troubleshooting Professional's dual nature. We've heard you, and starting next month you'll get two focused magazines, with the first Linux Productivity Magazine in August, and the next Troubleshooting Professional Magazine in September.
Steve
GNU/Linux is comprised of the Linux kernel originally crafted by Linus Torvalds, plus many, many utilities, a large number of which were utilities from the original GNU project. "GNU/Linux" is probably the most accurate moniker one can give to the operating system. Please be aware that in all of Troubleshooters.Com, when I say "Linux" I really mean "GNU/Linux". I completely believe that without the GNU project, without the GNU Manifesto and the GNU/GPL license it spawned, the operating system I abbreviate with the name "Linux" never would have happened.
You can read more about the GNU project and free software at the GNU website, a link to which is contained in this magazine's URL's section.
Also know that most of the time I use the terms "free software" and
"Open Source" interchangably. Although they are two separate movements
with very different motivations, licenses that are Open Source are usually
also free software and vice versa. The Open Source philosophy stresses
business practicality through source code availability and the free software
philosophy stresses freedom as an end in itself. Personally, I think both
are vitally important.
This all happened about a year ago. Soon after I created the proof of concept, starting with a simple dropped ball, and culminating in a moon and earth simulation which shows the moon to orbit the Earth in between 28 and 29 days, and then just for fun on day 60 an asteroid hits the moon and knocks both the moon and the earth, which is gravitationally attached to the earth, into a new trajectory.
So this issue of Troubleshooting Professional Magazine delves into simulation
of trajectories at a small enough distance to be considered flat. It can
be accomplished with any true C++, but this proof of concept uses the GNU
C++ that comes with any Linux distribution.
A vector is a measure of both magnetude and direction. Your car's velocity is a vector. If you're going north at 30MPH, the magnetude is 30MPH, and the direction is north. If you're going northeast at 64MPH, your magnetude is 64MPH, and your direction is northeast. The magnitude is often referred to as "the length" of the vector. Magnitudes typically represent speed, force, acceleration, or distance from another point.
Vectors are 3 dimensional (north/south, east/west, up/down), but when first learning about them it's easier to ignore one dimension and think of them in 2 dimensions.
64MPH northeast = 45.25MPH north + 45.25MPH east
Remembering high school geometry, this is nothing but an equilateral right triangle:
The preceding figure is an equilateral right triangle. It's a right triangle because one of the angles is a right angle (90 degrees). It's equilateral because each of the 2 sides adjacent to the right angle are equal. The hypotinuse (the line opposite the right ange) is 1.414 times the length of the sides adjacent to the right triangle. 1.414 is the square root of 2.
This is actually a special case of geometry for right trianges. If one calls the bottom side x, the top side y, and the hypotinuse h, then for any right triangle,
x2 + y2 = h2
Plugging in 1 for each of x and y,
12 + 12 = h2
So
h2 = 2
And therefore h is the square root of 2, or 1.414.
If the sides adjacent the right angle are 3 and 4, everything comes out even:
32 + 42 = 9 + 16 = 25 = h2
So h is the square root of 25, which is 5. This is called a "345" triangle:
Here's the general formula for the hypotenuse:
h = sqrt(x2 + y2)
where function sqrt() stands for "square root".
Likewise, solving for x yields the following:
x = sqrt(h2 - y2)
30MPH north * -1 = 30MPH south.
Obviously, dividing a vector by a scalar yields the same vector as multiplying the vector by the reciprocal of the scalar.
You add two or more vectors by summing their x values to make the new
vector's x, and the y values to make the new vector's y. This is typically
done when you add forces. Subtracting vector 1 from vector 2 is the same
as adding the inverse of vector 1 to vector 2. The inverse of vector 1
is a vector of the same magnitude but in the opposite direction.
| Multiplication | Vector | Scalar |
| Vector | N/A | The vector's magnitude is multiplied by the scalar |
| Scalar | The vector's magnitude is multiplied by the scalar | The product of the two numbers |
| Addition | Vector | Scalar |
| Vector | The result of adding two vectors is a second vector whose x is the
sum of the x values for both vectors, and whose y value is the sum of the
y values for both vectors. Can be represented graphically by placing the
addend vectors head to tail:
 |
N/A |
| Scalar | N/A | The sum of the two numbers |
| Subtraction | Vector | Scalar |
| Vector | Second vector whose x is the difference of the x values for both vectors,
and whose y value is the difference of the y values for both vectors. Can
be represented graphically by placing the tail of a directional inverse
of the vector being subtracted against the head of the other vector:
|
N/A |
| Scalar | N/A | The sum of the two numbers |
The trig functions work with right triangles. A right triangle has 3
sides and 3 angles. By definition of right triangle, one of the angles
is
90 degrees. The trig functions enable you to calculate any side or angle
based on any combination of 3 angles or sides, counting the right triangle.
The trig functions are simple ratios of various right triangle sides. These
ratios change depending on the angle under consideration. In the preceding
drawing, if we consider the angle in the lower left, the hypotinuse is
the blue line, the side adjacent is the line whose length is 4, and the
side opposite is the side whose length is 3. The following chart gives
definitions of the 6 trig functions. Note that the shaded functions are
less common, and can be easily deduced from the others:
| Name | Function symbol | Definition | |
| Sine | sin(angle) | Side opposite over hypotenuse | |
| Cosine | cos(angle) | Side adjacent over hypotinuse | |
| Tangent | tan(angle) | Side opposite over side adjacent | |
| Cotangent | cot(angle) | Side adjacent over hypotinuse | |
| Secant | sec(angle) | Hypotinuse over side adjacent (1/cosine) | |
| Cosecant | csc(angle) | Hypotinuse over side opposite (1/sine) |
There are tables giving these ratios for every angle. Better yet, standard C and C++ implement these functions in math.h.
h2 = x2 + y2 + z2
3 dimensions are hard to visualize on a computer monitor, which is why we've discussed the 2 dimensional oversimplification. But once you add the third dimension, everything's pretty straightforward.
////////////////////////////////////////////////////////
// vector.h, 3d vector class
#ifndef VECTOR_H
#define VECTOR_H
#include <math.h>
class Vector {
private:
double x,y,z;
public:
Vector(double xa,double ya,double za);
double getx() const {return x;};
double gety() const {return y;};
double getz() const {return z;};
static double getx(const Vector &v) {return v.x;};
static double gety(const Vector &v) {return v.y;};
static double getz(const Vector &v) {return v.z;};
void setx(const double xa) {x=xa;};
void sety(const double ya) {y=ya;};
void setz(const double za) {z=za;};
Vector vectorDistanceFrom(const Vector &v) const ;
double scalarDistanceFrom(const Vector &v) const ;
double getMagnitude() const {return(sqrt((x*x)+(y*y)+(z*z)));};
static double getMagnitude(const Vector &v) {
return(sqrt((v.x*v.x)+(v.y*v.y)+(v.z*v.z)));
};
static Vector addVectors(const Vector &v1,const Vector &v2) {
return(Vector(v1.getx()+v2.getx(),v1.gety()+v2.gety(),v1.getz()+v2.getz()));
}
void addVector(const Vector &v){
x+=v.getx();
y+=v.gety();
z+=v.getz();
}
Vector getUnitVector() const ;
static Vector getUnitVector(const Vector &v);
Vector timesScalar(const double d){x *= d; y *= d; z *= d;return(Vector(x,y,z));};
};
#endif
|
////////////////////////////////////////////////////////
// vector.cpp, 3d vector class
#include "vector.h"
Vector::Vector(double xa,double ya,double za) : x(xa), y(ya), z(za) {
}
Vector Vector::vectorDistanceFrom(const Vector &v) const {
return(Vector(v.x-x,v.y-y,v.z-z));
}
double Vector::scalarDistanceFrom(const Vector &v) const {
Vector vT=vectorDistanceFrom(v);
return(sqrt((vT.x * vT.x)+(vT.y * vT.y)+(vT.z * vT.z)));
}
Vector Vector::getUnitVector() const {
return(Vector(
x/getMagnitude(),
y/getMagnitude(),
z/getMagnitude()
));
}
Vector Vector::getUnitVector(const Vector &v) {
return(Vector(
v.x/getMagnitude(v),
v.y/getMagnitude(v),
v.z/getMagnitude(v)
));
}
|
Now let's test the class with the simplest possible test jig. This test
jig program adds two vectors with the same direction, and prints out the
magnitude and direction (as a unit vector). It does this by adding both
vectors into a total vector which is initially zero. Here's the test jig,
named vectortest.cpp:
////////////////////////////////////////////////////////
// vectortest.cpp, vector test jig
//
#include <stdio.h>
#include <iostream.h>
#include <math.h>
#include "vector.h"
void addstraight(void) {
Vector vec1(100,0,0);
Vector vec2(200,0,0);
Vector vecTotal(0,0,0);
vecTotal.addVector(vec1);
vecTotal.addVector(vec2);
cout << "Total magnitude is " << vecTotal.getMagnitude() << ".\n";
double totx = vecTotal.getUnitVector().getx();
double toty = vecTotal.getUnitVector().gety();
double totz = vecTotal.getUnitVector().getz();
cout << "Total direction x=" << totx << ", total direction y=" << toty << ", total direction z=" << totz << ".\n";
}
main() {
addstraight();
}
|
Here's the make file, called makefile.vec:
# make file for vectortest, a test jig for the Vector class vectortest: vectortest.o vector.o g++ -o vectortest vectortest.o vector.o; vector.o: vector.cpp body.h vector.h g++ -c -o vector.o vector.cpp; vectortest.o: vectortest.cpp vector.h g++ -c -o vectortest.o vectortest.cpp; clean: rm -rf vectortest.o vector.o |
Here's what happens when you run it:
[slitt@mydesk sim]$ make -f makefile.vec clean rm -rf vectortest.o vector.o [slitt@mydesk sim]$ make -f makefile.vec g++ -c -o vectortest.o vectortest.cpp; g++ -c -o vector.o vector.cpp; g++ -o vectortest vectortest.o vector.o; [slitt@mydesk sim]$ ./vectortest Total magnitude is 300. Total direction x=1, total direction y=0, total direction z=0. [slitt@mydesk sim]$ |
As expected, you get a magnitude that's the sum of the magnitudes of the original vectors, and a direction that's identical to that of the original vectors.
Now let's add an all-x vector magnitude 3, and an all-y vector magnitude
4, and see if we get a magnitude 5 vector in the proper direction, whose
x to y ratio is 3/4, or 0.75. Rather than writing the whole program, we'll
simply add a subroutine called add345, and change the subroutine call in
the main program:
void add345(void) {
Vector vec1(3,0,0);
Vector vec2(0,4,0);
Vector vecTotal(0,0,0);
vecTotal.addVector(vec1);
vecTotal.addVector(vec2);
cout << "Total magnitude is " << vecTotal.getMagnitude() << ".\n";
double totx = vecTotal.getUnitVector().getx();
double toty = vecTotal.getUnitVector().gety();
double totz = vecTotal.getUnitVector().getz();
cout << "Total direction x=" << totx << ", total direction y=" << toty << ", total direction z=" << totz << ".\n";
cout << "Direction X/Y is " << (float) totx/toty << ".\n";
}
|
Here's what happens when you run it:
[slitt@mydesk sim]$ ./vectortest Total magnitude is 5. Total direction x=0.6, total direction y=0.8, total direction z=0. Direction X/Y is 0.75. [slitt@mydesk sim]$ |
So much for 2 dimensional vectors. Let's add equal vectors in the x,
y and z directions, and see if the direction is equal in all 3 dimensons,
and if the resulting magnitude is sqrt(3) times the magnitude of each individual
vector. We make a subroutine called add3dim() to add an all-x, an all-y
and an all-z vector, and see what we get:
void add3dim(void) {
Vector vec1(1,0,0);
Vector vec2(0,1,0);
Vector vec3(0,0,1);
Vector vecTotal(0,0,0);
cout << "We\'re expecting an answer whose magnitude is ";
cout << "root of 3 (" << sqrt(3) << "),\n";
cout << "and whose direction is equal in all 3 directions.\n\n";
vecTotal.addVector(vec1);
vecTotal.addVector(vec2);
vecTotal.addVector(vec3);
cout << "Total magnitude is " << vecTotal.getMagnitude() << ".\n";
double totx = vecTotal.getUnitVector().getx();
double toty = vecTotal.getUnitVector().gety();
double totz = vecTotal.getUnitVector().getz();
cout << "Total direction x=" << totx << ", total direction y=" << toty << ", total direction z=" << totz << ".\n";
}
|
Here's what happens when you run it:
[slitt@mydesk sim]$ ./vectortest We're expecting an answer whose magnitude is root of 3 (1.73205), and whose direction is equal in all 3 directions. Total magnitude is 1.73205. Total direction x=0.57735, total direction y=0.57735, total direction z=0.57735. [slitt@mydesk sim]$ |
In the preceding example the total magnitude comes out to the square root of 3, which is exactly what you'd expect. The unit vector of this vector is 0.577 in each direction.
Finally, we'll make a Vector(2,3,6), which has a magnitude of 7 (22
+ 32 + 62 = 72), and then we'll multiply
it by 2:
void mult3dim(void) {
Vector vec1(2,0,0);
Vector vec2(0,3,0);
Vector vec3(0,0,6);
Vector vecTotal(0,0,0);
vecTotal.addVector(vec1);
vecTotal.addVector(vec2);
vecTotal.addVector(vec3);
cout << "Total magnitude is " << vecTotal.getMagnitude() << ".\n";
double totx = vecTotal.getUnitVector().getx();
double toty = vecTotal.getUnitVector().gety();
double totz = vecTotal.getUnitVector().getz();
cout << "Total direction x=" << totx << ", total direction y=" << toty << ", total direction z=" << totz << ".\n";
cout << "Vector components are x=" << vecTotal.getx() << ", ";
cout << "y=" << vecTotal.gety() << ", ";
cout << "z=" << vecTotal.getz() << ".\n";
cout << "\n";
cout << "Now multiply by 2:\n";
vecTotal.timesScalar(2);
cout << "Total magnitude is " << vecTotal.getMagnitude() << ".\n";
cout << "Total direction x=" << totx << ", total direction y=" << toty << ", total direction z=" << totz << ".\n";
cout << "Vector components are x=" << vecTotal.getx() << ", ";
cout << "y=" << vecTotal.gety() << ", ";
cout << "z=" << vecTotal.getz() << ".\n";
}
|
And here's how it looks when you run it:
[slitt@mydesk sim]$ make -f makefile.vec g++ -c -o vectortest.o vectortest.cpp; g++ -o vectortest vectortest.o vector.o; [slitt@mydesk sim]$ ./vectortest Total magnitude is 7. Total direction x=0.285714, total direction y=0.428571, total direction z=0.857143. Vector components are x=2, y=3, z=6. Now multiply by 2: Total magnitude is 14. Total direction x=0.285714, total direction y=0.428571, total direction z=0.857143. Vector components are x=4, y=6, z=12. [slitt@mydesk sim]$ |
If 2 bodies d Kilometers apart (at their centers) have masses m1 and m2 Kilograms respectively, the gravitational force between them is defined by this equation:
F = Gm1m2/d2
G is a constant 6.67e-11 with units meters2/killograms2.
For reasons beyond the scope of this article, the attraction is the same whether all the mass is compressed into the size of a marble, or whether it's as big as the earth, as long as there is some empty space between the two bodies. In other words, if the sun collapsed into a black hole the size of a pinhead, Newtonian physics says the earth would still maintain the same orbit as before. Relativity may change this, and also if the earth creates tides on the sun that would change in a solar collapse, but to a great degree the orbit would remain the same. Note that this entire tutorial is based only on Newtonian physics.
The most fundamental physics equation is:
F=ma
Where F is the total force on a body, m is the body's mass, and a is the body's acceleration. In gravitational problems, a body's mass cancels out of the equation. For instance, we all know gravitational force of a 1000Kg car is much more than a 1Kg basketball. That's why we can throw a basketball half way down the court, but we can't even lift the front end of the car. In fact, mass is proportional to the gravitational force. But it's also part of the right side of the equation, so the two can cancel out. In fact, all objects at the surface of the earth fall at an accelleration of 9.8 m/sec2.
Looking at it from a universe point of view, from ball's perspective:
Fball = Gmearthmball/d2
and Fball = mballaball
So substituting for force,
mballaball = Gmearthmball/d2
Dividing both sides by mball, we obtain the final acceleration equation:
aball = Gmearth/d2
Substituting G=6.67e-11, mearth=5.98e24, and d=6.37e6, (between the center of earth and the softball, i.e. the radius of the earth), we come up with this:
aball=6.67e-11 x 5.98e24/(6.37e6)2
= (6.67 x 5.98)e13/(6.37 x 6.37)e12
= 6.67 x 5.98 x 10/6.37 x 6.37
= 398.866/40.5769
= 9.8298
So everything's consistent beyond rounding errors. These body formulas are the basis of trajectory simulation. In the next article we'll discuss trajectory simulation.
But when you add in factors like air resistance or other bodies or forces, these equasions become significantly inaccurate unless they're modifed to become much more complex. Trajectory simulation, on the other hand, remains somewhat understandable even with these other factors. The technique is basically to timeslice the trajectory into tiny time increments. At each increment, you look at each body, sum all the forces upon it, and calculate its acceleration. From the acceleration you can deduce the new speed, and from the speed you can deduce the change in position. Yes, it's that simple.
Simulation expert Darren Humphrey from Simulation.Com suggested that I first calculate gravity effects on all bodies, then move all the bodies as required. I decided to do both in one step, and boy was I sorry. Doing each as a single step means that if there are a large number of bodies (say 1000), body 1 will encounter a very different environment than body 1000. It's not only inaccurate, but it's hard to code a single "move" step instead of a gravity calculation and then a movement.
So the correct (Darren's) process of trajectory simulation looks something like this:
timeIncrement=0.01
until sought event
foreach body
calculate effects from
all forces
foreach body
update acceleration
and velocity
move the body (velocity
times timeIncrement)
check for sought event
Of course, because we're representing phyical bodies, the only efficient way to code such a simulation is with objects. So in the second to next article we'll make a Body class with Vector variables to represent location and speed, and also has properties mass, diameter, name, color (in case you ever implement it graphically), and the universal gravitational constant G (6.67e-11 with units meters2/killograms2).
So the second to next article details the Body class. But before describing the Body class, we need to describe various helper classes to facilitate the Body class's functionality.
////////////////////////////////////////////////////////
// universeinfo.h, class to represent universal properties
#ifndef UNIVERSEINFO_H
#define UNIVERSEINFO_H
class UniverseInfo {
public:
const double G() const {return 6.67e-11;};
};
#endif
|
This class is instantiated as a global variable. Because it has no properties, its routines could also be called as class methods.
////////////////////////////////////////////////////////
// earthinfo.h, earth information class
#ifndef EARTHINFO_H
#define EARTHINFO_H
// moon info from http://www.solarviews.com/eng/moon.htm
class EarthInfo {
public:
double earthMass(){return 5.98e24;}; //kilograms
double G(){return 6.67e-11;};
double earthRadius(){return 6.37e6;}; //meters
double earthSurfaceAcceleration(){ return 9.8;}; //meters/sec^2
double moonMass(){return 7.349e22;}; //kilometers
double moonRadius(){return 1737400;}; //meters
double moonMeanOrbitalSpeed(){return 1030;}; //meters per second
double moonEquitorialSurfaceGravity(){return 1.62;}; //meters/sec^2
double moonMeanDistanceFromEarth(){return 3.844e8;}; //meters
double moonOrbitalPeriodInDays(){return 27.32166;}; //days
double moonRotationalPeriodInDays(){return 27.32166;}; //days
};
#endif
|
This class is instantiated as a global variable. Because it has no properties, its routines could also be called as class methods.
|
// conversions.h, conversion class #ifndef CONVERSIONS_H
#include <math.h> class Conversions{
#endif |
This class is instantiated as a global variable. Because it has no properties, its routines could also be called as class methods.
////////////////////////////////////////////////////////
// body.h, class to represent physical body
#ifndef BODY_H
#define BODY_H
#include "vector.h"
#include "universeinfo.h"
class Body {
Vector location;
Vector velocity;
double mass;
double diameter;
const char * name;
int color;
UniverseInfo uib;
Vector forceOnBody; //accumulator for external forces
//******* pfWindForce is a hook for including wind resistance ******
Vector (*pfWindForce) ( //callback calculating wind force
const Vector &, //body location
const Vector &, //body velocity
const double, //body's surface area
const double //"stickiness" of body's skin
);
public:
Body(const Vector &,const Vector &, const double, const double,const char *, const int);
Vector getLocation() const {return location;};
Vector getVelocity(){return velocity;};
const double getMass() const {return mass;};
const double getDiameter() const {return diameter;};
const char * getName() const {return name;};
const int getColor() const{return color;};
const double getx() const{return location.getx();};
const double gety() const{return location.gety();};
const double getz() const{return location.getz();};
const double getVelocityX() const{return velocity.getx();};
const double getVelocityY() const{return velocity.gety();};
const double getVelocityZ() const{return velocity.getz();};
void setLocation(const Vector &v) {location=v;};
void setVelocity(const Vector &v) {velocity=v;};
void setMass(const double m) {mass=m;};
Vector vectorDistanceFrom(const Body &b2) const { return getLocation().vectorDistanceFrom(b2.getLocation()); }
double scalarDistanceFrom(const Body &b2) const { return getLocation().scalarDistanceFrom(b2.getLocation()); }
void addGravityEffect(const Body&);
void updatePosition(double);
//*********** Specific to wind resistance **********
void addWindEffect();
void setWindForceCallback(Vector (*pf)(const Vector &, const Vector &, const double, const double)) {
pfWindForce=pf;
};
};
#endif
|
You'll notice most of the implementation is actually done in this .h file. Only the more complex functions, addGravityEffect(), updatePosition(), addWindEffect(), and the constructor are left to be defined in body.cpp.
The basic properties are self explanatory. location and velocity are vectors. mass and diameter are double precision floating point numbers. forceOnBody is a vector representation of the sum of all forces on the body at an instant in time, and is used to determine its movement (f=ma). The name and color aren't used in these exercises, but they become crucial if you hook this up to a realtime graphics system. The class provides all the set and get functions necessary to manipulate these properties in a proof of concept.
This class is remarkably robust in its support of wind resistance. This paragraph and the next 2 describe that support, but if it confuses you, you can skip it now and come back to it when needed.
A callback function pointer, called pfWindForce, can be called by the class (via the addWindEffect() method), to determine the wind-caused force on the body. This callback function is called with arguments corresponding to many of the properties that could have an effect in wind resistance, namely body location, body velocity, body's surface area and the "stickiness" of body's skin. Note that on anything but a spherical body, the surface area would be dependent on orientation, and that would need to be passed. But for the purpose of these exercises, we assume the bodies are spherical. In fact, we don't bother with surface area or skin stickiness, but instead use a constant available to the callback function that incorporates those properties. The body's velocity is always vital to the calculation of wind resistance. The location is important only if the trajectory goes through areas of varying air density. A ballistic missile would be such an application.
The addWindEffect() method does nothing unless a callback function has been installed with the setWindForceCallback() function. If so installed, addWindEffect() calls the callback to get the wind-caused force, and then adds it to the forceOnBody vector in much the same way that addGravityEffect() adds gravitational effects of other bodies.
The following is body.cpp, which implements some of the more
complex functions of the Body class. Remember that uib.G()
returns the universal gravitational constant:
////////////////////////////////////////////////////////
// body.cpp, class to represent physical body
#include "body.h"
#include <stdio.h>
Body::Body(
const Vector & initialLocation,
const Vector & initialVelocity,
const double massA,
const double diameterA,
const char * nameA,
const int colorA)
:
location(initialLocation),
velocity(initialVelocity),
mass(massA),
diameter(diameterA),
name(nameA),
color(colorA),
forceOnBody(Vector(0.0,0.0,0.0)),
pfWindForce(NULL){
};
void Body::addGravityEffect(const Body & otherBody) {
double distanceScalar = scalarDistanceFrom(otherBody);
double forceScalar = (uib.G()*mass*otherBody.getMass()/(distanceScalar * distanceScalar));
Vector forceVector = vectorDistanceFrom(otherBody);
forceVector=forceVector.getUnitVector();
forceVector.timesScalar(forceScalar);
forceOnBody.addVector(forceVector);
};
void Body::updatePosition(double timeInterval) {
// Calculate new acceleration
Vector accelVector=forceOnBody;
accelVector.timesScalar(1/mass);
double accelScalar=Vector::getMagnitude(accelVector);
// Calculate new velocity
Vector dVelocity=accelVector;
dVelocity.timesScalar(timeInterval);
Vector nextVelocity=Vector::addVectors(velocity,dVelocity);
// Calculate change in position
Vector medianVelocity=Vector::addVectors(velocity,nextVelocity);
medianVelocity.timesScalar(0.5);
Vector dLocation=medianVelocity;
dLocation.timesScalar(timeInterval);
// Update object properties
location.addVector(dLocation);
velocity=nextVelocity;
forceOnBody=Vector(0,0,0);
}
//*********** Specific to wind resistance **********
void Body::addWindEffect() {
if(pfWindForce) {
Vector windForceVector = pfWindForce(location,velocity,0,0);
forceOnBody.addVector(windForceVector);
}
}
|
In the preceding, the constructor (Body) simply does initialization. addGravityEffect() calculates the gravitational attraction to its Body argument, and adds that attractive force to the forceOnBody force accumulator. The addWindEffect() performs a similar function for wind forces, but only if a callback function has been defined to calculate wind force.
After all gravitational forces have been summed, and if appropriate the wind resistance force has been summed in, the net force on the object is known. Note that calculating the net force requires a call to addGravityEffect() for every other significant body in the system, and these forces must be totalled before a movement is calculated.
After the net force on all bodies in the system are known, those bodies must be moved accordingly. Each body is moved with a call to its updatePosition() method. This method first calculates an accelleration based on the net force, then calculates the new velocity, then calculates the body's new position, and finally updates the body's location and velocity properties. updatePosition() appears more complex than it is because vector math is done with unary functions, instead of arithmetic operators. C++ gurus can easily adapt operators (+, - and *) to handle vector math, after which pdatePosition() will appear trivial.
Note that addWindEffect() contains code preventing it from operating without an installed callback routine.
Now that you've defined the Body class, you need a test jig to exercise and prove it. The test jig is simply a ball being dropped from 100 meters, to see how long it will take to hit the earth. If you were to solve this by equation instead of by simulation, the equation is:
y = 1/2 * at2
In the preceding equation, y is the height dropped from, a is the acceleration, which should be about 9.8 meters/second2, and t is the time in seconds. Solve for T by first dividing each side by 1/2 * a, and then taking the square root of each side:
2y/a = t2
sqrt(2y/a) = t
Plugging in values 100 meters for y and 9.8 meters/sec2 for a,
t = sqrt(2*100/9.8) = sqrt(20.41) = 4.52 seconds
And since velocity is the initial velocity plus the product of acceleration and time, the velocity on impact is 4.52*9.8, or 44.27 meters/second2.
Now let's simulate it and see how close we come:
////////////////////////////////////////////////////////
// dropball.cpp, dropped ball class test jig
//
#include <stdio.h>
#include "body.h"
#include "vector.h"
#include "universeinfo.h"
#include "earthinfo.h"
EarthInfo ei;
UniverseInfo ui;
void showProgress(const Body & b1, double time) {
printf("%6.3f %6.3f %8.3f %8.3f\n",time,b1.getx(),b1.gety(), b1.getVelocityY());
}
void dropBall(){
const double heightDroppedInMeters = 100;
const double interval=0.25; //seconds between simulations
Body bBall(
Vector(0.0, heightDroppedInMeters, 0.0), // initial location
Vector(0.0, 0.0, 0.0), // velocity
0.5, // mass
0.1, // diameter
"ball", // name
0xff0000 // color (red)
);
Body bEarth(
Vector(0.0, -ei.earthRadius(), 0.0), // initial location
Vector(0.0, 0.0, 0.0), // velocity
ei.earthMass(), // mass
2.0 * ei.earthRadius(), // diameter
"earth", // name
0x00ff00 // color (green)
);
printf(" DROPPING BALL...\n");
printf(" TIME X Y SPEED\n");
double ellapsed=0.0;
showProgress(bBall, ellapsed);
while(bBall.gety() >= 0.0) {
bBall.addGravityEffect(bEarth);
bEarth.addGravityEffect(bBall);
bBall.updatePosition(interval);
bEarth.updatePosition(interval);
ellapsed += interval;
showProgress(bBall, ellapsed);
}
}
main() {
dropBall();
}
|
Follow the action in the dropBall subroutine. The height and time intervals are established, after which Body objects are instantiated for both the earth and the ball.
Next we begin our loop, which will become ever more familiar as you do the rest of the exercises in this Troubleshooting Professional issue. Each body's addGravityEffect() is called for the other body to add that gravitational attraction to the overall force on the body. Note that in this simple case there are only two bodies and no wind, so there is only one call to addGravityEffect() for each body. If you were evaluating an entire solar system there would be many more calls.
After the forces are added in, each body's updatePosition() method is called to update its position and velocity.
You compile it like this:
make dropballWhen it runs it will look something like this:
[slitt@mydesk sim]$ ./dropball DROPPING BALL... TIME X Y SPEED 0.000 0.000 100.000 0.000 0.250 0.000 99.693 -2.457 0.500 0.000 98.771 -4.915 0.750 0.000 97.235 -7.372 1.000 0.000 95.085 -9.830 1.250 0.000 92.321 -12.287 1.500 0.000 88.942 -14.744 1.750 0.000 84.948 -17.202 2.000 0.000 80.341 -19.659 2.250 0.000 75.119 -22.117 2.500 0.000 69.283 -24.574 2.750 0.000 62.832 -27.031 3.000 0.000 55.767 -29.489 3.250 0.000 48.087 -31.946 3.500 0.000 39.794 -34.404 3.750 0.000 30.886 -36.861 4.000 0.000 21.363 -39.319 4.250 0.000 11.226 -41.776 4.500 0.000 0.475 -44.233 4.750 0.000 -10.890 -46.691 [slitt@mydesk sim]$ |
The simulation shows the fall to take a little more than 4.5 seconds, and the finished velocity to be a little more than 44.23 meters per second. These figures are remarkably close to the figures predicted by the simple physics equations earlier in this article.
The concept has been proved, and you have simulated.
The following program uses the an earth Body object and a ball Body
object to simulate a throw in a world without air resistance. It will be
explained following the code.
////////////////////////////////////////////////////////
// throwball.cpp, ball throw, no wind test jig
//
#define MAIN
#include <stdio.h>
#include "body.h"
#include "vector.h"
#include "conversions.h"
#include "universeinfo.h"
#include "earthinfo.h"
EarthInfo ei;
UniverseInfo ui;
void showProgress(const Body & b1, double time) {
printf("%16.3f %16.3f %16.3f %16.3f\n",time,b1.getx(),b1.gety(), b1.getz());
}
void throwBall(){
const double interval=0.1; //seconds between simulations
const double throwSpeedMilesPerHour = 65;
const double throwAngleDegrees = 40.0;
const double throwHeightMeters = 1.5; //Height of cocked arm
double metersPerSecond =
Conversions::milesPerHourToMetersPerSecond(throwSpeedMilesPerHour);
Vector throwVelocity(0.0, cos(Conversions::degreesToRadians(throwAngleDegrees)),
sin(Conversions::degreesToRadians(throwAngleDegrees)));
throwVelocity.timesScalar(metersPerSecond);
printf("\nx is sidedrift, y is distance, z is height\n");
printf("throwVelocity=(%f,%f,%f)\n",throwVelocity.getx(),throwVelocity.gety(),throwVelocity.getz());
printf("\n Time SideDrift Distance Height\n\n");
Body bBall(
Vector(0.0, 0.0, throwHeightMeters), // initial location
throwVelocity, // velocity
0.5, // mass
0.1, // diameter
"ball", // name
0xff0000 // color (red)
);
Body bEarth(
Vector(0.0,0.0,-ei.earthRadius()), // initial location
Vector(0.0,0.0,0.0), // velocity
ei.earthMass(), // mass
2.0 * ei.earthRadius(), // diameter
"earth", // name
0x00ff00 // color (green)
);
double ellapsed=0.0;
showProgress(bBall, ellapsed);
while(bBall.getz() >= 0.0) {
bBall.addGravityEffect(bEarth);
bEarth.addGravityEffect(bBall);
bBall.updatePosition(interval);
bEarth.updatePosition(interval);
ellapsed += interval;
showProgress(bBall, ellapsed);
}
}
main() {
throwBall();
}
|
The simulation logic is contained in throwBall(). The simulation sampling interval, the throw velocity, the throw angle, and the height of the thrower's hand are established in the correct units. The throw velocity is a vector. In this exercise the y axis is assumed to be the direction of the throw, the z axis is vertical, and the x axis is horizontal and perpendicular to the throw.
Next a 1/4 Kg ball is instantiated as a Body object, and the Earth is instantiated as a Body object.
Now we go into the familiar loop where the ball and earth calculate the net force on themselves as a result of each other, and then both positions are updated, the ellapsed time is updated, and the new positions are shown. The loop terminates when z goes below 0. Z is height, and 0 height is the earth's surface, so it terminates when the ball hits the ground.
Notice how easy this was. It's almost identical to the dropped ball, except that the ball starts in a different place with a different velocity. There's a place to plug in almost any fact about the bodies being simulated.
You compile this program with this command:
make throwballThe following is the output that results from running the program:
x is sidedrift, y is distance, z is height throwVelocity=(0.000000,22.254655,18.673873) Time SideDrift Distance Height 0.000 0.000 0.000 1.500 0.100 0.000 2.225 3.318 0.200 0.000 4.451 5.038 0.300 0.000 6.676 6.660 0.400 0.000 8.902 8.183 0.500 0.000 11.127 9.608 0.600 0.000 13.353 10.935 0.700 0.000 15.578 12.163 0.800 0.000 17.804 13.294 0.900 0.000 20.029 14.325 1.000 0.000 22.255 15.259 1.100 0.000 24.480 16.094 1.200 0.000 26.706 16.831 1.300 0.000 28.931 17.470 1.400 0.000 31.157 18.010 1.500 0.000 33.382 18.452 1.600 0.000 35.607 18.796 1.700 0.000 37.833 19.041 1.800 0.000 40.058 19.189 1.900 0.000 42.284 19.237 2.000 0.000 44.509 19.188 2.100 0.000 46.735 19.040 2.200 0.000 48.960 18.794 2.300 0.000 51.186 18.450 2.400 0.000 53.411 18.007 2.500 0.000 55.637 17.466 2.600 0.000 57.862 16.827 2.700 0.000 60.087 16.090 2.800 0.000 62.313 15.254 2.900 0.000 64.538 14.320 3.000 0.000 66.764 13.287 3.100 0.000 68.989 12.157 3.200 0.000 71.215 10.928 3.300 0.000 73.440 9.600 3.400 0.000 75.666 8.175 3.500 0.000 77.891 6.651 3.600 0.000 80.117 5.029 3.700 0.000 82.342 3.308 3.800 0.000 84.567 1.489 3.900 0.000 86.793 -0.428 |
This follows the same trajectory curve you learned in college physics. You'll note that a 65 MPH throw (fast, but not professional athlete speed), results in an extremely long 86.7 meter throw (282 feet). This is how far you would throw in a vacuum sealed football field. As you'll see in the next exercise, that number goes way down in the air.
Before leaving this exercise, try modifying the throw angle and speed. You'll find that in a vacuum, the optimum throwing angle is the 45 degrees predicted in high school physics. This too will change when we work with air resistance.
Conceptually adding air resistance and wind effect is simple. With each iteration you simply call the ball's addWindEffect() right after its addGravityEffect() in order to add wind resistance to the total force load on the ball. The trick is calculating the wind effect. That calculation is done by a callback function completely separated from the Body object, and connected to the Body object through a pointer at runtime. That function is passed everything about the Body object that is relevant to air resistance. In this very simple simulation it's passed the body's location, velocity, surface area and surface stickiness.
The function is charged with using these, and any other info it might have, to calculate the wind resistance. In our case we will create a Windinfo class and a global Windinfo object for this function to work with and store info in.
Now that you know the theory, these are the steps to convert the airless throw simulation to a wind-aware throw simulation:
This ball throw is perhaps the simplest possible wind-aware simulation. It's reasonable to assume that the wind will be uniform throughout the travel of the throw, and throughout the time the ball is in the air. It's not like a cannon shot where the wind velocity changes every 300 meters or so, or like a rocket reentry where air density varies. Because the ball is a sphere, its orientation is not an issue, unlike a real life cannon shot. It's so simple that we ignore the "stickiness" and surface area, and simply assign a constant that yields reasonable results. That constant is 0.0015.
But it's still a little tricky. One might for a minute suppose that
on a windless day these wind calculations would not come into play. But
in fact every moment of the ball's travel, its own velocity relative to
the still air creates a retarding force. On a windy day, the wind velocity
is vector subtracted from the ball velocity to calculate the headwind.
That headwind is then squared and multiplied by the constant to get the
force. The direction of the force is the direction of the headwind.
| ////////////////////////////////////////////////////////
// windinfo.h, wind class and callback for ball throw #ifndef WINDINFO_H
class WindInfo {
#ifdef MAIN
Vector cbWindForce(
//callback calculating wind force
//******* Calculate headwind ******
#endif
|
The Windinfo object simply stores the wind velocity and the constant to be used in wind force calculations. This global object provides a way you can pass info into the callback function. There are probably more OOP ways of doing this, but I was pressed for time, it works, and it's modular. So when the callback function needs the wind velocity, it looks at wi.getWindVelocity(). And when it needs the constant, it gets wi.getWindCoefficient(). Because this object is global, the main routine can change these values.
Start with the cbWindForce() function in the preceding code. This functions return type is that expected by the Body class, and its argument types are those expected by the Body class. When a Body object calls it, the Body object passes this function the body's location, velocity, surface area and stickiness. The latter two are not used by this function.
The function subtracts any wind velocity from the body's velocity to calculate the headwind. It then extracts the headwind's magnitude, squares it, multiplies it by the constant, and then multiplies that value into a unit vector derived from the headwind. The result is a force proportional to the square of the headwind, which most people agree is a reasonable way to predict wind force.
Notice the code that actually instantiates the WindInfo object:
#ifdef MAIN
WindInfo wi(Vector(0.0, 0.0, 0.0));
#endif
As long as MAIN is defined in the main program, that code will
put the global variable in the main file, and only the main file, which
is what you want. So object wi will be available and in scope
in the main file.
|
Was this callback function a good idea? It requires a global variable (the Windinfo object) to communicate between the program and the callback function. Its arguments are cast in concrete by the coding of the Body class, and these arguments are woefully inadequate for many complex wind calculations. It probably would have been better to make a class to represent wind force calculations, subclass it as necessary (with late binding, of course), and give the Body class a pointer to it. Or perhaps subclass the Body class itself to cover various situations. After all, orientation, profile and the like are traits of the Body, not of the wind. The bottom line is I had a limited time to create this tutorial, and the callback function is very modular and adaptable, so I used it. |
bBall.setWindForceCallback(cbWindForce); //windforce wi.setWindCoefficient(0.0015); //light, but reasonable //windforce wi.setWindVelocity(Vector(1,-0.5,0).getUnitVector().timesScalar(Conversions::milesPerHourToMetersPerSecond(0))); //windforce |
bBall.addWindEffect(); //windforce |
#2 first associates the callback function in windinfo.h, cbWindForce(), as the Body object bBall's callback function. Then the second line sets the wind constant to a reasonable value. The third line sets the wind velocity. The third line is set up so you define the wind speed in the argument to milesPerHourToMetersPerSecond(), and you define the direction in the Vector() initialization. The reason the Vector() defines only a direction is because its unit vector is extracted.
#3 simply tells bBall to add the wind resistance to other forces acting on it.
The following is the complete source code for throwball_wind.cpp:
////////////////////////////////////////////////////////
// throwball_wind.cpp, ball throw, wind included test jig
//
#define MAIN
#include <stdio.h>
#include "body.h"
#include "vector.h"
#include "conversions.h"
#include "universeinfo.h"
#include "earthinfo.h"
#include "windinfo.h" //windforce
EarthInfo ei;
UniverseInfo ui;
void showProgress(const Body & b1, double time) {
printf("%16.3f %16.3f %16.3f %16.3f\n",time,b1.getx(),b1.gety(), b1.getz());
}
void throwBall(){
const double interval=0.1; //seconds between simulations
const double throwSpeedMilesPerHour = 65;
const double throwAngleDegrees = 40.0;
const double throwHeightMeters = 1.5; //Height of cocked arm
double metersPerSecond =
Conversions::milesPerHourToMetersPerSecond(throwSpeedMilesPerHour);
Vector throwVelocity(0.0, cos(Conversions::degreesToRadians(throwAngleDegrees)),
sin(Conversions::degreesToRadians(throwAngleDegrees)));
throwVelocity.timesScalar(metersPerSecond);
printf("\nx is sidedrift, y is distance, z is height\n");
printf("throwVelocity=(%f,%f,%f)\n",throwVelocity.getx(),throwVelocity.gety(),throwVelocity.getz());
printf("\n Time SideDrift Distance Height\n\n");
Body bBall(
Vector(0.0, 0.0, throwHeightMeters), // initial location
throwVelocity, // velocity
0.25, // mass
0.1, // diameter
"ball", // name
0xff0000 // color (red)
);
bBall.setWindForceCallback(cbWindForce); //windforce
wi.setWindCoefficient(0.0015); //light, but reasonable //windforce
wi.setWindVelocity(
Vector(1,-0.5,0).getUnitVector().timesScalar(Conversions::milesPerHourToMetersPerSecond(0))
);
Body bEarth(
Vector(0.0,0.0,-ei.earthRadius()), // initial location
Vector(0.0,0.0,0.0), // velocity
ei.earthMass(), // mass
2.0 * ei.earthRadius(), // diameter
"earth", // name
0x00ff00 // color (green)
);
double ellapsed=0.0;
showProgress(bBall, ellapsed);
while(bBall.getz() >= 0.0) {
bBall.addGravityEffect(bEarth);
bBall.addWindEffect(); //windforce
bEarth.addGravityEffect(bBall);
bBall.updatePosition(interval);
bEarth.updatePosition(interval);
ellapsed += interval;
showProgress(bBall, ellapsed);
}
}
main() {
throwBall();
}
|
The preceding had the same 65mph throw velocity, 40 degree release angle, and 1.5 meter throw height as the no-air example in the preceding article. The wind speed is 0 (wi.setWindVelocity()). But the air has some viscosity (wi.setWindCoefficient()), and that retards the throw significantly. Now the throw goes only 63 meters instead of the 85 meters when thrown in a vacuum. Remember, the wind speed has been set to 0.
The following is the output of the program:
x is sidedrift, y is distance, z is height throwVelocity=(0.000000,22.254655,18.673873) Time SideDrift Distance Height 0.000 0.000 0.000 1.500 0.100 0.000 2.206 3.302 0.200 0.000 4.374 4.975 0.300 0.000 6.507 6.522 0.400 0.000 8.606 7.946 0.500 0.000 10.673 9.250 0.600 0.000 12.710 10.437 0.700 0.000 14.718 11.509 0.800 0.000 16.698 12.468 0.900 0.000 18.653 13.315 1.000 0.000 20.582 14.054 1.100 0.000 22.487 14.685 1.200 0.000 24.370 15.211 1.300 0.000 26.230 15.631 1.400 0.000 28.069 15.949 1.500 0.000 29.888 16.165 1.600 0.000 31.686 16.280 1.700 0.000 33.465 16.296 1.800 0.000 35.225 16.213 1.900 0.000 36.967 16.032 2.000 0.000 38.690 15.756 2.100 0.000 40.395 15.384 2.200 0.000 42.082 14.917 2.300 0.000 43.751 14.357 2.400 0.000 45.403 13.705 2.500 0.000 47.037 12.961 2.600 0.000 48.654 12.127 2.700 0.000 50.253 11.204 2.800 0.000 51.834 10.193 2.900 0.000 53.398 9.095 3.000 0.000 54.943 7.911 3.100 0.000 56.471 6.643 3.200 0.000 57.980 5.292 3.300 0.000 59.471 3.858 3.400 0.000 60.943 2.345 3.500 0.000 62.397 0.752 3.600 0.000 63.832 -0.919 |
#5 is interesting because the wind badly overpowers the throw. The throw goes only 15.9 meters, and during the last .3 seconds of the throw the ball actually travels back toward the thrower. Experiment with the throw angle, and note that just like in real life, you do better throwing lower into a headwind. In fact, with this particular throw the optimal angle of release is roughly 18 degrees.
You might also try playing with the release angle with tailwinds, and notice that the optimal angle is a little more than 45 degrees, but very little more. As a matter of fact, you'll notice that in most normal conditions, a release angle of 40 degrees gives a nice, long throw.
Try varying the mass of the ball, and notice that heavier balls do better in severe headwinds, just as you would expect.
If you really want to have fun, change the program so throw velocity, release angle, wind velocity, and wind constant are arguments to throwBall(), and use a loop to vary them. For instance, you could use a double loop to vary headwind, and within that release angle, coming up with a list of optimal release angles to use at every windspeed.the loop to vary release angle, thereby
We've got one more simulation, but it's a doosey, and if it were written accurately, it would also be somewhat useful on a battlefield. But first, some theory...
This magazine features proofs of concept, so absolute accuracy isn't needed. In situations like real battlefield or strategic applications of artillery and intercontinental ballistic missiles, accuracy must be kept to the ultimate achievable. This rules out much approximation. As a matter of fact, it rules out all significant approximation when calculating the landing point based on a launch point, launch direction, launch elevation, and wind conditions. All approximations will cause error in the final landing point.
But there's still a role for approximations in simulation.
So far we've dealt with problems where the launch velocity vector and wind conditions were known, and the target was being solved for. Many actual problems have a known target and wind conditions, and the launch velocity vector is being solved for. This is done by a trial and error (negative feedback) algorithm like the following:
|
.------------------.
| | V | ---------------- | | Take a shot | | ---------------- | | | V | ---------------------- | | Evaluate the error | | ---------------------- | | | V | --------------------- | | Correct the | | | launch velocity | | --------------------- | | | `------------------' |
Any approximation in taking the shot will cause errors in the launch velocity calculations. But the next 2 steps, evaluating the error and making the correction, are merely feedback, and as such they can contain approximations as long as:
#2 is necessary so that you don't oscillate around the correct answer, but instead steadily home in on it. If this proves difficult with the math knowledge at hand, one can always brute force it by substituting more tiny change increments for a good calculation, and you'll very slowly but very surely converge.
If this all sounds like I'm using simulation as a crutch for inadequate math knowledge, you're absolutely right. But I'm in good company. Ultimately most systems could be predicted by equations, but doing so is beyond current human knowledge. So we humans have developed ways to substitute simulation for equations.
This TPM issue makes many approximations. That's reasonable for a proof of concept. If you decide to create a program to aim cannons or missiles on the battlefield, you'd hire a physics PH.D from UCF or MIT to give you the equations yielding the best approximations, thereby reducing the number of trials to a minimum.
A proof of concept is a simple implementation whose purpose is learning. This article details a proof of concept cannon aimer that works through a negative feedback mechanism whereby successive iterations correct for errors in previous iterations.
As a proof of concept, this is not only a simplification, but in fact it has bugs. When you finish the exercises in this article you will note that the example, as given, fails when the target is exactly due East of the cannon. I didn't debug further because I ran out of time.
So as you use this proof of concept, learn from it, find its flaws, improve it, and if you're really interested, try improving its basic design. You would not use this cannon aimer on the battlefield. But as a proof of concept -- a learning tool, it performs quite well. Once you understand the code and the principles behind it, you might want to redesign the planShot_low() and planShot_high() functions, because when code performs like unruly hair usually indicates a bad design. I leave it as an exercise to the astute reader to design this program to work consistently without all sorts of exception code.
Whether or not you choose to accept this challenge, if you carefully read this article and work with its code, trajectory simulation will never again seem like magic.
This cannon aiming computer kept doing successive simulations until it simulated a direct hit, and then delivered the direction to aim, and how high to aim.
The top of the program looks like this:
//////////////////////////////////////////////////////// // cannon.cpp, find cannon vector muzzle angle test jig // #define MAIN #include <stdio.h> #include <stdlib.h> #include <string.h> #include <math.h> #include "body.h" #include "vector.h" #include "conversions.h" #include "universeinfo.h" #include "earthinfo.h" #include "windinfo.h" //windforce EarthInfo ei; // Earth measurements and constants UniverseInfo ui; // Universal measurements and constants |
In the preceding, we define MAIN for program-global variables. stdio.h, stdlib.h, string.h and math.h are needed for input/output, process, string manipulation and math functions such as trig, logs and exponents.
The need for body.h, vector.h, conversions.h, universeinfo.h, earthinfo.h, and windinfo.h were discussed in the ball throw with wind resistance simulation.
To return the result of a test shot, we define a struct Snapshot
like this:
struct Snapshot {
double t;
double x;
double y;
double z;
};
|
A struct Snapshot is a simple way to show the x, y and z location of the vector, plus the flight time.
A class called Givens serves up the given information in the
problem of cannon aiming. It looks like this:
// *** Given information expressed as a global object
class Givens {
public:
const double maxRangeInMiles(){return(2);};
const double rangeAngleInDegrees(){return(40);};
const double rangeWindspeed(){return(0);};
const double windCoefficient(){return(.001);};
const double shellDiameterInInches(){return(8);};
const double shellWeightInPounds(){return(100);};
const double shellMass(){return(shellWeightInPounds() * Conversions::kilogramsPerPound());};
const double maxRange(){
return(Conversions::metersPerKilometer()* Conversions::kilometersPerMile() * maxRangeInMiles());
};
const double rangeAngleInRadians(){return(rangeAngleInDegrees()/Conversions::degreesPerRadian());};
const double shellDiameter(){return(shellDiameterInInches()/(12.0 * Conversions::feetPerMeter()));}
};
Givens givens;
|
As you can see in the preceding code, a global instance of Givens is instantiated as givens, so that givens can be queried by any code in this file.
For simulating a cannon shot, we'll make a function called shootCannon()
that's very similar to the ball throwing programs discussed previously
in this TPM issue. Here's the shootCannon() function:
// *** Based on launch point, vector muzzle velocity, target height,
// *** and simulation interval, plot the shot.
// *** showTrajectory is a flag indicating whether to print the
// *** entire trajectory, instead of just the hit point.
const struct Snapshot shootCannon(
Vector launchPoint,
Vector muzzleVelocity,
const double targetHeight,
const double interval,
int showTrajectory
){
if(showTrajectory) {
printf("%16s %16s %16s %16s\n",
"Time",
"X",
"Y",
"Z");
}
Body bMissile(
launchPoint, // initial location
muzzleVelocity, // velocity
givens.shellMass(), // mass
givens.shellDiameter(), // diameter
"shell", // name
0xff0000 // color (red)
);
bMissile.setWindForceCallback(cbWindForce); //windforce
wi.setWindCoefficient(givens.windCoefficient());
Body bEarth(
Vector(0.0,0.0,-ei.earthRadius()), // initial location
Vector(0.0,0.0,0.0), // velocity
ei.earthMass(), // mass
2.0 * ei.earthRadius(), // diameter
"earth", // name
0x00ff00 // color (green)
);
double ellapsed=0.0;
int crossedOver = 0;
if(targetHeight <= bMissile.getz())
crossedOver++; //start from above target height
Vector prevPosition(0,0,0);
double prevTime=0;
while(crossedOver < 2) {
// Update break logic
prevTime = ellapsed;
prevPosition=Vector(bMissile.getx(), bMissile.gety(), bMissile.getz());
// Calculate forces
bMissile.addGravityEffect(bEarth);
bMissile.addWindEffect(); //windforce
bEarth.addGravityEffect(bMissile);
// Update velocities and positions based on forces
bMissile.updatePosition(interval);
bEarth.updatePosition(interval);
ellapsed += interval;
if(showTrajectory) {
printf("%16.3f %16.3f %16.3f %16.3f\n",
ellapsed,
bMissile.getx(),
bMissile.gety(),
bMissile.getz());
}
// Detect over to under and under to over transitions.
if(crossedOver == 0) {
if(bMissile.getz() > targetHeight) {
crossedOver++;
}
else if (prevPosition.getz() > bMissile.getz()) { // already falling
crossedOver = 101; // never reaches target height
}
}
else {
if(bMissile.getz() < targetHeight) {
crossedOver++;
}
}
}
// Return the hit point and flight time
struct Snapshot s;
if(crossedOver > 100) { //never reached height of target
s.x = launchPoint.getx();
s.y = launchPoint.gety();
s.z = launchPoint.getz();
s.t = 0;
}
else {
// Interpolate based on height above and below
double tooFarFraction = (targetHeight-bMissile.getz())/(bMissile.getz() - prevPosition.getz());
s.x = bMissile.getx() + (bMissile.getx() - prevPosition.getx()) * tooFarFraction;
s.y = bMissile.gety() + (bMissile.gety() - prevPosition.gety()) * tooFarFraction;
s.z = bMissile.getz() + (bMissile.getz() - prevPosition.getz()) * tooFarFraction;
s.t = ellapsed - (ellapsed - prevTime) * tooFarFraction;
}
return(s);
}
|
The preceding code is so similar to the ball throwing routines described
earlier in this magazine that all we need discuss is its arguments:
| Argument | Usage |
| launchPoint | This is the Vector location of the cannon. |
| muzzleVelocity | This is the vector velocity of the shell as it emerges from the cannon. |
| targetHeight | The height at which you consider the shell to have hit. Obviously, the shell strikes higher points before it strikes lower ones. For practical purposes, this is used to aim at targets above or below the cannon. |
| interval | The time intervals at which forces are analyzed and locations and velocities recalculated. Smaller figures yield more accurate results, but require greater computing power. |
| showTrajectory | A flag. If set, the shell's progress is listed at each interval period. |
The return value is a struct Snapshot, which contains the hit point's X, Y and Z value, as well as the shell's flight time.
Now that we have a way of predicting the landing point of a shell shot with a particular velocity, the next step is to find the muzzle speed which yields a 2 mile shot at 40 degrees on a windless day. We'll call the subroutine getMuzzleSpeed(), and it will return the scalar muzzle speed. That subroutine is basically a looped ball throw designed to show us the muzzle speed that will go 2 miles at a 40 degree elevation angle on a windless day. We'll use that muzzle velocity for the aiming function.
When I say basically, we modify the algorithm to accommodate a performance
consideration. With any iterative approximation, there's a tradeoff between
small iterations and computing time. This is unfortunate, because smaller
iterations yield much more accurate results because linear interpolation
is inaccurate with non-linear functions such as a falling body. We get
the best of both worlds by getting an approximate answer with large iterations,
then using a little bit less than that number as a starting point and using
smaller iterations, on and on until further shrinking of the difference
between iterations falls below a threshold of acceptability:
double getMuzzleSpeedRough(double low, double interval) {
// Find transition from higher to lower than ground
double angle=givens.rangeAngleInRadians();
struct Snapshot s;
s.y=0;
struct Snapshot sPrev;
double speed=low; // meters per second
double prevSpeed= 0;
while(s.y < givens.maxRange()) {
prevSpeed = speed;
sPrev.t=s.t; sPrev.x=s.x; sPrev.y=s.y; sPrev.z=s.z;
s=shootCannon(Vector(0,0,0), Vector(0, speed*cos(angle), speed*sin(angle)), 0, 0.1, 0);
speed += interval;
}
// Interpolate to exact z=0
double tooFarFraction = (s.y - givens.maxRange())/(s.y - sPrev.y);
double finalSpeed = speed - (speed - prevSpeed) * tooFarFraction;
// Return
return(finalSpeed);
}
double getMuzzleSpeed() {
double speed = 0;
double interval=10;
double prevSpeed=-10;
while((speed-prevSpeed)*(speed-prevSpeed) > .0000000001) {
prevSpeed = speed;
speed = getMuzzleSpeedRough(speed - 12*interval, interval);
interval /= 10;
}
return(speed);
}
|
In the preceding, getMuzzleSpeedRough() starts at argument low with interval interval, and repeatedly does shootCannon() at increasing angles until the projectile goes farther than 2 miles (givens.maxRange()). At that time there's an interpolation between the last shot to fall short and the first shot to fall long, and the interpolated speed is passed back.
getMuzzleSpeed() repeatedly calls getMuzzleSpeedRough() with more accurate starting points and smaller intervals, until two successive getMuzzleSpeedRough() calls produce results within .00001 meters of each other (the square of the difference is .0000000001).
You might wonder why we subtract 12*interval from the starting speed in the getMuzzleSpeedRough() call. A casual reading of the program indicates that 2 would be a sufficient multiple, but in fact 10 is required for the program not to diverge in the long direction. A closer look at the program reveals that in the call to getMuzzleSpeedRough(), the interval value is 1/10 of the value of interval in the previous iteration, and the previous iteration's value of interval reveals the maximum divergence from the desired value. Therefore you must multiply interval by 10. I used 12 just to provide an added margin of error in case of some odd non-linear effects. Using 12 doesn't add much compute time.
Over and over again throughout the cannon aimer program you'll see similar algorithms, where there's a rough calculator that repeatedly gets called by a looping calculator continually passing ever finer intervals and more accurate starting points until two successive calls fall within the required level of accuracy.
We also need a function to verify the accuracy of getMuzzleSpeed().
That function is called
// *** Confirms muzzle speed produces a 2 mile shot at 40 degrees elevation on a windless day
void confirmMuzzleSpeed(const double muzzleSpeed) {
printf("y should be %f miles, or %f meters\n", givens.maxRangeInMiles(), givens.maxRange());
double angle = givens.rangeAngleInRadians();
struct Snapshot s=shootCannon(
Vector(0,0,0),
Vector(0, muzzleSpeed*cos(angle),
muzzleSpeed*sin(angle)),
0,
0.1,
0
);
printf("Shot Confirmation: speed=%8.2f, t=%8.2f, x=%8.2f, y=%12.6f, z=%12.6f\n",
muzzleSpeed, s.t, s.x, s.y, s.z);
}
|
To facilitate human understanding of directions, and limit the likelihood of getting arguments wrong, the program includes a direction to degrees-from-due-east converter called findAngle(). This function gives us the polar coordinate angle represented by an English description of an angle. For example,
double myAngle = findAngle(20, "degrees", "west", "of", "north")
printf("myAngle=%f\n", myAngle);
The preceding code would print 110, which is the 90 degrees to straight
north, plus the 20 degrees west of that. The following is the code for
findAngle(),
plus the code for findAngleUsage(), which is a syntax printing
routine called when wrong arguments are passed to
findAngle():
// *** Syntax description for findAngle() function
void findAngleUsage() {
fprintf(stderr, "findAngle usage:\n");
fprintf(stderr, "findAngle(degrees, \"degrees\", direction, \"of\", direction)\n");
fprintf(stderr, "possible combinations are:\n");
fprintf(stderr, "\"north\" \"of\" \"east\" \n");
fprintf(stderr, "\"east\" \"of\" \"north\" \n");
fprintf(stderr, "\"west\" \"of\" \"north\" \n");
fprintf(stderr, "\"north\" \"of\" \"west\" \n");
fprintf(stderr, "\"south\" \"of\" \"west\" \n");
fprintf(stderr, "\"west\" \"of\" \"south\" \n");
fprintf(stderr, "\"east\" \"of\" \"south\" \n");
fprintf(stderr, "\"south\" \"of\" \"east\" \n");
fprintf(stderr, "\n\n");
exit(1);
}
// *** Converts an English description of a direction into a polar coordinate angle
double findAngle(
double degrees,
const char * degreesT, // Must be the literal "degrees"
const char * direction1,
const char * ofT, // Must be the literal "of"
const char * direction2
) {
if(strcasecmp(degreesT, "degrees")) findAngleUsage();
if(strcasecmp(ofT, "of")) findAngleUsage();
double angle = 0;
if ((!strcasecmp(direction1, "north"))&&(!strcasecmp(direction2, "east"))) {angle = degrees;}
else if((!strcasecmp(direction1, "east"))&&(!strcasecmp(direction2, "north"))) {angle = 90-degrees;}
else if((!strcasecmp(direction1, "west"))&&(!strcasecmp(direction2, "north"))) {angle = 90+degrees;}
else if((!strcasecmp(direction1, "north"))&&(!strcasecmp(direction2, "west"))) {angle = 180-degrees;}
else if((!strcasecmp(direction1, "south"))&&(!strcasecmp(direction2, "west"))) {angle = 180+degrees;}
else if((!strcasecmp(direction1, "west"))&&(!strcasecmp(direction2, "south"))) {angle = 270-degrees;}
else if((!strcasecmp(direction1, "east"))&&(!strcasecmp(direction2, "south"))) {angle = 270+degrees;}
else if((!strcasecmp(direction1, "south"))&&(!strcasecmp(direction2, "east"))) {angle = 360-degrees;}
else {findAngleUsage();}
return(angle);
}
|
The program also includes a test jig function to test the findAngle() function. The test jig is called testAngles():
// *** Test jig for the findAngle() function
void testAngles(void) {
printf("%5.2f\n", findAngle(15, "degrees", "north", "of", "east"));
printf("%5.2f\n", findAngle(15, "degrees", "east", "of", "north"));
printf("%5.2f\n", findAngle(15, "degrees", "west", "of", "north"));
printf("%5.2f\n", findAngle(15, "degrees", "north", "of", "west"));
printf("%5.2f\n", findAngle(15, "degrees", "south", "of", "west"));
printf("%5.2f\n", findAngle(15, "degrees", "west", "of", "south"));
printf("%5.2f\n", findAngle(15, "degrees", "east", "of", "south"));
printf("%5.2f\n", findAngle(15, "degrees", "south", "of", "east"));
printf("%5.2f\n", findAngle(15, "degrees", "frunk", "of", "east"));
}
|
Next we'll load the WindInfo object with all its data, including the coefficient. It would be a couple hour project to change the WindInfo object to accept wind velocities at different elevations and interpolate between them. Doing so would make the program much more useful on the battlefield, but we'll forgo the realism in order to make a simple proof of concept.
With any shot short of maximum range there are two trajectories that can hit a target. One is low and relatively straight, and the other is high and comes in from above. The loop will be capable of calculating either, depending on whether it starts with a high elevation angle or a low one.
With any iterative approximation like this, there's a tradeoff between small iterations and computing time. This is unfortunate, because smaller iterations yield much more accurate results because linear interpolation is inaccurate with non-linear functions such as a falling body. In the getMuzzleSpeedRough() code we used a loop that repeatedly called a rough loop. That won't work here, because errors can be either on the long-short axis, or on the left-right axis. It would get way too complicated using the 2 subroutine method, with its repeated calculation of how much to lop off the angle. Fortunately there's another way to converge in while reducing the interval value. It's kind of slick. But first, a word on high and low trajectories.
There are two ways you can find the high trajectory:
As the target goes closer to the outer edge of the cannon's range, the
low and high trajectories get closer together. On any point exactly on
the cannon's range (after factoring in wind), the two trajectories converge
into 1. The following picture shows 2 trajectories much closer together
(meaning the target is much closer to the edge of the cannon's range):
With a severe headwind, both low and high trajectories will be below
45 degrees. With a huge tailwind, the low trajectory will be very low,
while the high trajectory is very high. Given enough of a tailwind, the
high trajectory could actually be more than 90 degrees, meaning you're
shooting away from the target but the wind blows it back to the target.
As you can imagine, the high trajectory is much more susceptable to crosswinds than the low trajectory. For that reason the high trajectory algorithm is much more likely to fail than the low trajectory algorithm. When would you use the high trajectory? It's very useful when shooting over high walls or hills. The low trajectory might crash into the hill, but the high trajectory flies over the hill and hits the target behind the hill.
For left-right aiming, we simply calculate the angle of divergence between target and hit point, and compensate by an amount half that angle. In that way we quickly converge on the target in a left-right context. Actually, in many cases we can divide by a number less than 2, and it will converge quickly. But if the divisor gets too small, it can oscillate around the correct value and send the program into an infinite loop. I found 1.65 to be a quick divisor that didn't seem to cause severe oscillation, but I chose to use 2 to be on the safe side. Using numbers like 3 or 4 just make the aiming process slower and sometimes cause a "timeout".
Note that for each iteration of short-long aiming, we complete a loop of one or more left-right aimings so at all times the hit point is close to the line between cannon and target. This prevents situations, at high wind velocity and ranges close to the cannon's maximum, in which the left-right angle is always playing catch-up with the elevation angle, and when it catches up it sends the hit point farther from the target. This results in oscillation and therefore many more iterations and slower aiming. In extreme cases it can cause the elevation angle continuously higher until the aimer terminates on a 180 degree angle, or worst case it goes past the maximum iterations.
So basically, left-right aiming iteratively closes the divergence angle between the target and the hit point.
Short-long isn't quite so simple. We have no easy way of translating the error into an angle, so we must simply increase elevation angle at a specific interval. When the hit crosses from too short to too long, we reverse the interval direction and divide it by 3. Then, when it crosses back to too-short, reverse and divide by 3 again. In this way both error and interval quickly shrink as accuracy becomes tight. In hindsight, this would have been a much better algorithm for getMuzzleSpeed() also, but we'll leave that routine as-is just to demonstrate there's more than one way to skin a cat.
Here's simplified pseudo-code of the aimer, with all angles described
in degrees:
zAngle = -89 xyAngle = targetAngle interval = 8 //number of degrees to raise cannon each try while(closeness > 1cm) prevOvershoot = overshoot while(abs(xyAngle-targetAngle) > interval/10) mV=findVectorMuzzleVelocity(muzzleSpeed, zAngle, xyAngle) hitPoint=shootCannon(mV, targetElevation) xyAngle -= (hitPoint.xyAngle - targetAngle)/2 closeness = scalar(hitPoint - target) overshoot = componentTowardTarget(hitPoint) - scalar(target) if(prevOvershoot * overshoot < 0) interval /= -3 zAngle += interval if((zAngle > 180) || (zAngle < -180)) break if(iteration > maxIterations) break |
// *** This is the cannon aimer for lower trajectories.
// *** It works by repeatedly calculating the error, re-aiming accordingly, and taking another shot
const Vector planShot_low(
const Vector target, // vector location of the target
const double muzzleSpeed // scalar speed of the cannon shell
) {
// Declare and init vars
Vector hitPoint(0,0,0); // Where the shell crosses the target height going down
Vector prevHitPoint(0,0,0); // Hit point on previous iteration (break logic)
double xyAngleHitPoint; // xy angular direction of hitpoint
double xyAngleDeviation; // xy angular deviation between target and hitpoint
Vector muzzleVelocity(0,0,0); // Vector muzzle velocity producing the current hit point
double zAngle = -89; // elevation angle, start very low for shooting into valley
double xyAngle; // directional angle
double closeness = 100000000; // scalar distance between hit point and target,
// start huge to start loop
double overshoot = -1; // difference between shot distance and target distance,
// start with undershoot
double prevOvershoot; // previous iteration's overshoot
double interval = 8; // amount to increase elevation angle each iteration
double intervalFactor = 3.0; // divisor to decrease interval on new overshoots or undershoots
int maxIterations = 2000; // point to assume you're infinitely looping
int iteration; // iteration number
int blownBack; // flag -- was it blown by high wind backward away from target
// Find angular location of target
double xyAngleTarget = atan2(target.gety(), target.getx()) * Conversions::degreesPerRadian();
xyAngle = xyAngleTarget;
// Repeatedly estimate and shoot until you hit within 1 cm (.01) of target
iteration=0;
while(closeness > .01) {
// set previous values to current
prevHitPoint = Vector(hitPoint.getx(), hitPoint.gety(), hitPoint.getz());
prevOvershoot = overshoot;
// For this elevation angle, make sure the xy angular deviation is less/eq the interval/10
xyAngleDeviation = 99999;
while(xyAngleDeviation > fabs(interval)/10.0) {
if(iteration++ > maxIterations) break;
// printf("diag %d: zAngle=%10.6f, xyAngle=%10.6f, blBk=%d, closeness=%9.5f\n",
// iteration, zAngle, xyAngle, blownBack, closeness);
// Translate angles to a velocity vector
double zt = muzzleSpeed * sin(zAngle/Conversions::degreesPerRadian());
double xyt = muzzleSpeed * cos(zAngle/Conversions::degreesPerRadian());
double xt = xyt * cos(xyAngle/Conversions::degreesPerRadian());
double yt = xyt * sin(xyAngle/Conversions::degreesPerRadian());
muzzleVelocity = Vector(xt, yt, zt);
// Take the shot
struct Snapshot s = shootCannon(Vector(0,0,0), muzzleVelocity, target.getz(), 0.1, 0);
hitPoint = Vector(s.x, s.y, s.z);
// Calculate the xy angle of deviation between hitpoint and target
xyAngleDeviation = 0;
double xyAngleCorrection = 0;
xyAngleHitPoint = xyAngleTarget;
blownBack = 0;
// don't correct if shot never achieved altitude
if(hitPoint.getMagnitude() > 0.00000000001) {
xyAngleHitPoint = atan2(hitPoint.gety(), hitPoint.getx()) * Conversions::degreesPerRadian();
xyAngleDeviation = xyAngleHitPoint - xyAngleTarget;
xyAngleCorrection = xyAngleDeviation/2;
if((xyAngleDeviation > 90) || (xyAngleDeviation < -90)) blownBack = 1;
}
// Adjust xy angle for next iteration
if(!blownBack) xyAngle -= xyAngleCorrection;
}
//Below this point work with elevation angle incrementation and loop termination
// Determine how close you got (scalar)
closeness = target.scalarDistanceFrom(hitPoint);
// Calculate how much short or long
double tempAngle = (xyAngleHitPoint - xyAngleTarget)/Conversions::degreesPerRadian();
overshoot = hitPoint.getMagnitude() * cos(tempAngle) - target.getMagnitude();
// Reverse and shorten interval on over to undershoot transition or vice versa
if(prevOvershoot * overshoot < 0) {
interval /= intervalFactor;
interval *= -1;
}
// Increment zAngle
zAngle += interval;
// The rest of the loop detects and acts on loop termination criteria
// Detect total inability to achieve the target at any trajectory
if((zAngle > 180) || (zAngle < -180)) {
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
// Detect and brake out of infinite loop
if(iteration > maxIterations) {
fprintf(stderr, "Program error: Hung while trying to calculate shot for\n");
fprintf(stderr, "target (%8.3f, %8.3f, %8.3f).\n",
target.getx(), target.gety(), target.getz());
fprintf(stderr, "Call programmer: internal error\n\n");
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
}
return(muzzleVelocity);
}
|
At the most basic level, the preceding code implements a nested loop, with the inner loop adjusting the left-right aiming (xyAngle), and the outer loop adjusting the muzzle elevation angle (zAngle). It's actually possible to use a single loop and adjust both in a single pass, but at the edge of the cannon's range at high winds the single loop algorithm can fail because xyAngle trails zAngle and falls hopelessly behind.
This loop is your basic negative feedback mechanism. First aiming angles xyAngle and zAngle are converted to a vector velocity called muzzleVelocity. Then a shot is simulated with that velocity. Finally, the left right deviation of the new shot from target is calculated, along with the new zAngle and xyAngle. Once the angular deviation from target is less than 1/10 of the zAngle interval, the inner loop terminates and passes control to the tail end of the outer loop.
One more thing. Did you notice a variable called blownBack?
The blownBack variable is a flag that is set if and only
if a headwind component manages to thrust the shell behind the cannon,
from the point of view of the target. See the pictures below:
In shot 1 above, a severe southwest wind blows the shell back south behind the initial firing line, and also causes it to drift significantly west. In shot 2 above, an almost completely southerly headwind blows the shell back behind the firing line. In each case, the angle from hitpoint to target is more than 90 degrees. When the angle is more than 90 degrees, you have blowback.
You don't want to recalculate xyAngle when there's blowback. Blowback is primarily a sign that your zAngle is too high for wind conditions. Recalculating xyAngle would cause excessive oscillation in the calculations, possibly leading to a complete divergence from the target. Blowback is handled by considering the shot an undershoot, leading to a decreased zAngle. Once blowback is eliminated by lowering zAngle, xyAngle can be adjusted.
You'll notice in the preceding code that if a blowback condition is reached in the inner loop, the inner loop will infinitely loop until rescued by exceeding maxIterations, at which time the calculation will fail. That's a bug. It's probably soluble simply by breaking the inner loop if blownBack is set, but I didn't have time to test all the possible side effects of such a change. The existing code works in a wide variety of target distances, target elevations, and wind conditions. It's relatively easy to code a cannon aimer that works in reasonable conditions, but when you test 300 meter/sec winds near the maximum cannon range, things get wonky.
There's a variable called closeness that represents the magnitude of the vector distance between the target and where the shell pierces a plane parallel to the ground at the altitude of the the target. For practical purposes, this is "how close the shell came" to the target. The loop test is closeness > .01. In other words, until the shot comes within 1 cm of the target, keep going. Once we drop through the inner loop, the first thing we do is calculate closeness.
The next thing we do is calculate whether it's an overshoot or an undershoot by determining whether the shot's component in the direction of the target is greater than the distance from cannon to target. On an overshoot condition, variable overshoot is a positive number. It's negative on undershoots. Blowbacks are undershoots, naturally.
If you remember, every time we toggle between overshoot and undershoot, we need to reverse the interval direction, and cut its magnitude by a factor of 3. The next few lines do that by multiplying overshoot by prevOvershoot. If the current and past outer loop iterations were undershoots, the product of two negatives is positive, and we do not cut and reverse. Likewise, if the current and past were overshoots, the product of two positives is positive, and we do not cut and reverse. But if the last iteration was an overshoot and this one is an undershoot, or if it's the other way around, the product of a positive and negative number is negative, causing the cutting and reversing of the interval variable.
Now that we know how much to increment, we increment zAngle in preparation for the next iteration. The final several lines of the outer loop detect termination conditions. Obviously the termination condition we'd like to see is the shell within a centimeter of the target, but sometimes that doesn't work out. If the target is out of range, zAngle will keep increasing until it goes past 180, in which case special code detects that condition and terminates the loop by setting closeness to 0, and signifying the error by setting muzzleVelocity to Vector(0,0,0). Sometimes, in very windy conditions and a target very close to the cannon's maximum range, the angle may go up for awhile, and then go down past -180, in which case it's terminated and the algorithm admits defeat.
The next few code lines bail if the number of inner iterations exceeds maxIterations. The maxIterations variable is sized at 2000 to allow all reasonable calculations, yet prevent infinite loops. Failures to converge routinely manifest themselves as infinite loops.
Another detail is the calculation of overshoot. This algorithm uses the component of the cannon to hitpoint distance on the cannon to target axis. I had originally used the scalar distance of the shot's hit point from the cannon, minus the scalar distance of target to cannon. They both worked, and seemed to work equally as well. I chose the component to prevent a case where an extreme crosswind would baloon the scalar shot distance, when in fact the cannon elevation could not have pushed the shell to the target.
Of course one could argue the other case. In a strong crosswind, after an 8 degree increase, the shot is blown right and thus the the component along the target vector is shorter than the target vector, but if the shot had been aimed left a degree more, it would have been right on the line, past the target. However, the double looping algorithm minimizes the likelihood of this happening.
All sorts of details crop up when the target is near the maximum range of the cannon. For instance, it's likely that a zAngle increment could go past the target, past the maximum distance, past the high trajectory, and therefore not overshoot. Under such a circumstance, the algorithm will keep incrimenting until either it goes past 180 degrees or it exceeds maxIterations, either of which will report an inability to reach the target. But in fact, if the intervals had been smaller, the target would have been found. This can be guarded against by a few lines of code which reverse and divide the interval if a higher shot produces a more negative overshoot. Problem is, a lag in the xyAngle could conceivably produce a loss of forward progress even below the zAngle of the low trajectory. Perhaps two consecutively more negative overshoots on two consecutively raised zAngle values would be better to trigger the reversal. Of course, if the target is right on the maximum range, the reverse interval might once again skip over the trajectory without converging. The bottom line is that this particular algorithm can't find the cannon aim necessary to hit a target right on the maximum range for a given wind condition. And in fact there's a finite band inside the maximum range in which this algorithm is likely to mistakenly conclude that the shot can't be made.
If you really want to succeed on targets very near the maximum range, perhaps you could use an algorithm that calculates the rate of change of overshoot against the rate of change of angle, on both sides of the maximum range, extrapolates both, and considers the intersection to be the maximum range. From there zAngle could be backed down slowly until the target is hit.
Indeed, a knowledge of mathematics far greater than mine would be a real asset in this cannon aiming simulation. If you really need to implment an accurate system, you'd be well advised to hire a Physics Ph.D from UCF or MIT to give your programmer ideas on how to best accommodate the many tricky situations of cannon aiming.
There's one more detail to consider. The cannon aimer is only as accurate as the shootCannon() function. If shootCannon() reports wrong hit points, the aimer will dutifully adjust muzzle velocity accordingly. Note the distinction. At all but the most extreme ranges and wind conditions, the cannon aimer is a negative feedback mechanism. By their very nature, negative feedback systems all but eliminate errors in their components. Thus you can use 2 or 3 or 4 as the xyAngle correction divisor, and the result will be equally accurate. But the shootCannon() routine is not a negative feedback mechanism, and is not included in the negative feedback of the aimer. It's simply an input. If this were a non-switching voltage regulator, shootCannon() would correspond to the zener diode on which all the output voltage is based.
So for ultimate accuracy, shootCannon() would need to be made more accurate. One way to do that would be to call it with ever decreasing time intervals until two consecutive results came back within some fraction of each other. Perhaps that fraction could be defined as the scalar difference between the two hit points, and perhaps when that distance falls to the same distance as the permissible miss by the aimer (1 cm in this example), the result would be considered accurate. Needless to say, this would increase the computing power required to aim the cannon. Luckily, heavy duty processors are getting cheaper all the time, and these algorithms could probably be adapted to Linux clusters.
I chose coming up from the low trajectory. The reason was simple. Vertical shots are at the mercy of the winds, and therefore pinning down a left-right correction on high shots is difficult. Keeping those corrections in sync with the descending vertical angle can be impossible in high crosswinds.
Contrast that with starting at the low trajectory. You're starting with a known good point, so you can use tiny intervals and move up, keeping directions and angles in sync at every step. Here's the pseudocode:
zAngle = zAngleLowTrajectory xyAngle = xyAngleLowTrajectory zAngle = addSlightAmount(zAngle) interval = 1 //number of degrees to raise cannon each try while(closeness > 1cm) prevOvershoot = overshoot while(abs(xyAngle-targetAngle) > interval/10) mV=findVectorMuzzleVelocity(muzzleSpeed, zAngle, xyAngle) hitPoint=shootCannon(mV, targetElevation) xyAngle -= (hitPoint.xyAngle - targetAngle)/2 closeness = scalar(hitPoint - target) overshoot = componentTowardTarget(hitPoint) - scalar(target) if(prevOvershoot * overshoot < 0) interval /= -3 zAngle += interval if((zAngle > 180) || (zAngle < -180)) break if(iteration > maxIterations) break |
You basically get far enough above the low trajectory elevation that closeness will be more than the 1cm necessary to terminate the loop, and then do the same thing as always -- toggle between overshoot and undershoot while narrowing the interval, all the while reigning in any left-right error. The one thing that's different here is that you shoot higher on overshoots, and lower on undershoots. Because the switch is a simple toggle, that's taken care of by the initial settings.
Because this code is so much like the low trajectory finder, there's no need of discussion other than to remember that high trajectory shots are much more susceptible to winds than their low trajectory brethren, which is why we set the initial interval to 1 degree instead of 8. For shots very near the cannon's range, you might want to change interval to 0.1, although on far range shots that might cause termination on maxIterations for shots that could have been made, so you might want to increase maxIterations and take the performance hit.
// *** This is the cannon aimer for higher trajectories.
// *** It works by repeatedly calculating the error, re-aiming accordingly, and taking another shot
// *** Arg muzzleVelocityA must be the muzzle velocity that produces the low trajectory hit.
// *** Therefore this function can only be called after a call to planShot_low().
const Vector planShot_high(
const Vector target, // vector location of the target
const Vector muzzleVelocityA // VECTOR speed of the cannon shell
) {
// Guarantee overshoot
struct Snapshot s;
Vector muzzleVelocity = muzzleVelocityA;
double muzzleSpeed = muzzleVelocity.getMagnitude();
double desiredCloseness = .01; // stop when you're this close
double closeness = 0; // scalar distance between hit point and target
double prevCloseness;
while(closeness <= desiredCloseness) {
muzzleVelocity.addVector(Vector(0,0,1));
s = shootCannon(Vector(0,0,0), muzzleVelocity, target.getz(), 0.1, 0);
closeness = target.scalarDistanceFrom(Vector(s.x, s.y, s.z));
}
double overshoot = 1; // difference between shot distance and target distance
// Calculate some values
Vector hitPoint = Vector(s.x, s.y, s.z);
double xyAngle = atan2(muzzleVelocity.gety(), muzzleVelocity.getx()) * Conversions::degreesPerRadian();
double xySpeedTemp = sqrt(muzzleVelocity.getx() * muzzleVelocity.getx() +
muzzleVelocity.gety() * muzzleVelocity.gety());
double zAngle = atan2(muzzleVelocity.getz(), xySpeedTemp) * Conversions::degreesPerRadian();
double xyAngleTarget = atan2(target.gety(), target.getx()) * Conversions::degreesPerRadian();
// Declare and init other vars
Vector prevHitPoint(s.x, s.y, s.z); // Hit point on previous iteration (break logic)
double prevOvershoot; // previous iteration's overshoot
double interval = 1; // amount to increase elevation angle each iteration
double intervalFactor = 3.0; // divisor to decrease interval on new over or undershoots
int maxIterations = 2000; // point to assume you're infinitely looping
int iteration; // iteration number
int blownBack = 0; // Blown behind launch.
double xyAngleDeviation; // xy angular deviation between target and hitpoint
double xyAngleHitPoint; // angle of the hit point (for xy angle corrections)
// Repeatedly estimate and shoot until you hit within 1 cm (.01) of target
iteration=0;
while(closeness > desiredCloseness) {
xyAngleDeviation = 999999;
while(xyAngleDeviation > fabs(interval)/10) {
if(iteration++ > maxIterations) break;
// printf("diag %4d: zAngle=%10.6f, xyAngle=%10.6f, closeness=%9.5f, blBk=%d\n",
// iteration, zAngle, xyAngle, closeness, blownBack);
// set previous values to current
prevHitPoint = Vector(hitPoint.getx(), hitPoint.gety(), hitPoint.getz());
prevOvershoot = overshoot;
prevCloseness = closeness;
// Take the shot
struct Snapshot s = shootCannon(Vector(0,0,0), muzzleVelocity, target.getz(), 0.1, 0);
hitPoint = Vector(s.x, s.y, s.z);
// Calculate the xy angle of deviation between hitpoint and target
xyAngleHitPoint = xyAngleTarget;
xyAngleDeviation = 0;
double xyAngleCorrection = 0;
blownBack = 0;
// Correct left-right if shot achieved altitude
if(hitPoint.getMagnitude() > -0.00000000001) {
xyAngleHitPoint = atan2(hitPoint.gety(), hitPoint.getx()) * Conversions::degreesPerRadian();
xyAngleDeviation = xyAngleHitPoint - xyAngleTarget;
xyAngleCorrection = xyAngleDeviation/2;
if((xyAngleDeviation > 90) || (xyAngleDeviation < -90)) blownBack = 1;
}
// Adjust xy angles for next iteration
if(!blownBack) xyAngle -= xyAngleCorrection;
// Translate angles to a new velocity vector
double zt = muzzleSpeed * sin(zAngle/Conversions::degreesPerRadian());
double xyt = muzzleSpeed * cos(zAngle/Conversions::degreesPerRadian());
double xt = xyt * cos(xyAngle/Conversions::degreesPerRadian());
double yt = xyt * sin(xyAngle/Conversions::degreesPerRadian());
muzzleVelocity = Vector(xt, yt, zt);
// Determine how close you got (scalar)
closeness = target.scalarDistanceFrom(hitPoint);
}
// Calculate how much short or long
if(blownBack) {
overshoot = -1; // If it's blown back, it's definitely an undershoot
}
else
{
double tempAngle = (xyAngleHitPoint - xyAngleTarget)/Conversions::degreesPerRadian();
overshoot = hitPoint.getMagnitude() * cos(tempAngle) - target.getMagnitude();
}
if(prevOvershoot * overshoot < 0) { // if switched between over and under or vice versa
interval /= intervalFactor;
interval *= -1;
}
// The rest of the loop detects and acts on loop termination criteria
zAngle += interval; // raise elevation angle
// Detect total inability to achieve the target at any trajectory
if((zAngle > 180) || (zAngle < -180)) {
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
// Detect and brake out of infinite loop
if(iteration > maxIterations) {
fprintf(stderr, "Program error: Hung while trying to calculate shot for\n");
fprintf(stderr, "target (%8.3f, %8.3f, %8.3f).\n",
target.getx(), target.gety(), target.getz());
fprintf(stderr, "Call programmer: internal error\n\n");
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
}
return(muzzleVelocity);
}
|
This routine outputs the vector muzzle velocity of both the low trajectory
and the high trajectory, and the travel time of each.
// *** This is a test jig for the planShot() function
void attempt(
const double muzzleSpeed, //scalar speed of shell on firing
const double degrees,
const char * degreesT, // literal "degrees"
const char * direction1,
const char * ofT, // literal "of"
const char * direction2,
const double targetHeight // height of target
) {
double targetAngle = findAngle(degrees, "degrees", direction1, "of", direction2);
printf("\n\nangle is %f\n\n", targetAngle);
Vector target(
cos(targetAngle/Conversions::degreesPerRadian()),
sin(targetAngle/Conversions::degreesPerRadian()),
1
);
target.timesScalar(1.4142135 * Conversions::kilometersPerMile()* Conversions::metersPerKilometer());
target.setz(targetHeight);
Vector muzzleVelocity = planShot_low(target, muzzleSpeed);
struct Snapshot s= shootCannon(Vector(0,0,0), muzzleVelocity, targetHeight, 0.1, 0);
double windDirection = atan2(wi.getWindVelocity().gety(), wi.getWindVelocity().getx());
windDirection *= Conversions::degreesPerRadian();
printf("\nWind velocity=%9.4f at %f degrees\n", wi.getWindVelocity().getMagnitude(), windDirection);
printf(" Proof: Low shot... t=%8.3f x=%8.3f y=%8.3f z=%8.3f\n",s.t,s.x,s.y, s.z);
printf("Muzzle velocity=(%7.2f, %7.2f, %7.2f), zAngle=%8.3f\n",
muzzleVelocity.getx(),
muzzleVelocity.gety(),
muzzleVelocity.getz(),
Conversions::degreesPerRadian() * atan2(muzzleVelocity.getz(),
sqrt(muzzleVelocity.getx()*muzzleVelocity.getx() + muzzleVelocity.gety()*muzzleVelocity.gety())));
printf("\n-----------------------------------------------------------------\n");
muzzleVelocity = planShot_high(target, muzzleVelocity);
s= shootCannon(Vector(0,0,0), muzzleVelocity, targetHeight, 0.1, 0);
printf("Proof High shot... t=%8.3f x=%8.3f y=%8.3f z=%8.3f\n",s.t,s.x,s.y, s.z);
printf("Muzzle velocity=(%7.2f, %7.2f, %7.2f), zAngle=%8.3f\n",
muzzleVelocity.getx(),
muzzleVelocity.gety(),
muzzleVelocity.getz(),
Conversions::degreesPerRadian() * atan2(muzzleVelocity.getz(),
sqrt(muzzleVelocity.getx()*muzzleVelocity.getx() + muzzleVelocity.gety()*muzzleVelocity.gety())));
}
|
Once you have attempt() coded, your main routine can look like this:
void main() {
wi.setWindVelocity(Vector(0,0,0));
double muzzleSpeed = getMuzzleSpeed();
confirmMuzzleSpeed(muzzleSpeed);
wi.setWindVelocity(Vector(14.14, -14.14, 0));
attempt(muzzleSpeed, 0, "degrees", "north", "of", "east", 777);
attempt(muzzleSpeed, 45, "degrees", "north", "of", "east", 777);
attempt(muzzleSpeed, 90, "degrees", "north", "of", "east", 777);
attempt(muzzleSpeed, 45, "degrees", "north", "of", "west", 777);
attempt(muzzleSpeed, 0, "degrees", "north", "of", "west", 777);
attempt(muzzleSpeed, 45, "degrees", "south", "of", "west", 777);
attempt(muzzleSpeed, 90, "degrees", "south", "of", "east", 777);
attempt(muzzleSpeed, 45, "degrees", "south", "of", "east", 777);
printf("\n\n");
}
|
The preceding first calculates the muzzle speed that can go 2 miles from a 40 degree trajectory on a windless day, then verifies the accuracy of that muzzle speed. Next it sets the wind velocity to 20 meters/second to the southeast (sqrt(14.142 + 14.142) is 20 meters per second.
////////////////////////////////////////////////////////
// cannon.cpp, find cannon vector muzzle angle test jig
//
#define MAIN
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "body.h"
#include "vector.h"
#include "conversions.h"
#include "universeinfo.h"
#include "earthinfo.h"
#include "windinfo.h" //windforce
EarthInfo ei; // Earth measurements and constants
UniverseInfo ui; // Universal measurements and constants
struct Snapshot {
double t;
double x;
double y;
double z;
};
// *** Given information expressed as a global object
class Givens {
public:
const double maxRangeInMiles(){return(2);};
const double rangeAngleInDegrees(){return(40);};
const double rangeWindspeed(){return(0);};
const double windCoefficient(){return(.001);};
const double shellDiameterInInches(){return(8);};
const double shellWeightInPounds(){return(100);};
const double shellMass(){return(shellWeightInPounds() * Conversions::kilogramsPerPound());};
const double maxRange(){
return(Conversions::metersPerKilometer()* Conversions::kilometersPerMile() * maxRangeInMiles());
};
const double rangeAngleInRadians(){return(rangeAngleInDegrees()/Conversions::degreesPerRadian());};
const double shellDiameter(){return(shellDiameterInInches()/(12.0 * Conversions::feetPerMeter()));}
};
Givens givens;
// *** Based on launch point, vector muzzle velocity, target height,
// *** and simulation interval, plot the shot.
// *** showTrajectory is a flag indicating whether to print the
// *** entire trajectory, instead of just the hit point.
const struct Snapshot shootCannon(
Vector launchPoint,
Vector muzzleVelocity,
const double targetHeight,
const double interval,
int showTrajectory
){
if(showTrajectory) {
printf("%16s %16s %16s %16s\n",
"Time",
"X",
"Y",
"Z");
}
Body bMissile(
launchPoint, // initial location
muzzleVelocity, // velocity
givens.shellMass(), // mass
givens.shellDiameter(), // diameter
"shell", // name
0xff0000 // color (red)
);
bMissile.setWindForceCallback(cbWindForce); //windforce
wi.setWindCoefficient(givens.windCoefficient());
Body bEarth(
Vector(0.0,0.0,-ei.earthRadius()), // initial location
Vector(0.0,0.0,0.0), // velocity
ei.earthMass(), // mass
2.0 * ei.earthRadius(), // diameter
"earth", // name
0x00ff00 // color (green)
);
double ellapsed=0.0;
int crossedOver = 0;
if(targetHeight <= bMissile.getz())
crossedOver++; //start from above target height
Vector prevPosition(0,0,0);
double prevTime=0;
while(crossedOver < 2) {
// Update break logic
prevTime = ellapsed;
prevPosition=Vector(bMissile.getx(), bMissile.gety(), bMissile.getz());
// Calculate forces
bMissile.addGravityEffect(bEarth);
bMissile.addWindEffect(); //windforce
bEarth.addGravityEffect(bMissile);
// Update velocities and positions based on forces
bMissile.updatePosition(interval);
bEarth.updatePosition(interval);
ellapsed += interval;
if(showTrajectory) {
printf("%16.3f %16.3f %16.3f %16.3f\n",
ellapsed,
bMissile.getx(),
bMissile.gety(),
bMissile.getz());
}
// Detect over to under and under to over transitions.
if(crossedOver == 0) {
if(bMissile.getz() > targetHeight) {
crossedOver++;
}
else if (prevPosition.getz() > bMissile.getz()) { // already falling
crossedOver = 101; // never reaches target height
}
}
else {
if(bMissile.getz() < targetHeight) {
crossedOver++;
}
}
}
// Return the hit point and flight time
struct Snapshot s;
if(crossedOver > 100) { //never reached height of target
s.x = launchPoint.getx();
s.y = launchPoint.gety();
s.z = launchPoint.getz();
s.t = 0;
}
else {
// Interpolate based on height above and below
double tooFarFraction = (targetHeight-bMissile.getz())/(bMissile.getz() - prevPosition.getz());
s.x = bMissile.getx() + (bMissile.getx() - prevPosition.getx()) * tooFarFraction;
s.y = bMissile.gety() + (bMissile.gety() - prevPosition.gety()) * tooFarFraction;
s.z = bMissile.getz() + (bMissile.getz() - prevPosition.getz()) * tooFarFraction;
s.t = ellapsed - (ellapsed - prevTime) * tooFarFraction;
}
return(s);
}
double getMuzzleSpeedRough(double low, double interval) {
// Find transition from higher to lower than ground
double angle=givens.rangeAngleInRadians();
struct Snapshot s;
s.y=0;
struct Snapshot sPrev;
double speed=low; // meters per second
double prevSpeed= 0;
while(s.y < givens.maxRange()) {
prevSpeed = speed;
sPrev.t=s.t; sPrev.x=s.x; sPrev.y=s.y; sPrev.z=s.z;
s=shootCannon(Vector(0,0,0), Vector(0, speed*cos(angle), speed*sin(angle)), 0, 0.1, 0);
speed += interval;
}
// Interpolate to exact z=0
double tooFarFraction = (s.y - givens.maxRange())/(s.y - sPrev.y);
double finalSpeed = speed - (speed - prevSpeed) * tooFarFraction;
// Return
return(finalSpeed);
}
double getMuzzleSpeed() {
double speed = 0;
double interval=10;
double prevSpeed=-10;
while((speed-prevSpeed)*(speed-prevSpeed) > .0000000001) {
prevSpeed = speed;
speed = getMuzzleSpeedRough(speed - 12*interval, interval);
interval /= 10;
}
return(speed);
}
// *** Confirms muzzle speed produces a 2 mile shot at 40 degrees elevation on a windless day
void confirmMuzzleSpeed(const double muzzleSpeed) {
printf("y should be %f miles, or %f meters\n", givens.maxRangeInMiles(), givens.maxRange());
double angle = givens.rangeAngleInRadians();
struct Snapshot s=shootCannon(
Vector(0,0,0),
Vector(0, muzzleSpeed*cos(angle),
muzzleSpeed*sin(angle)),
0,
0.1,
0
);
printf("Shot Confirmation: speed=%8.2f, t=%8.2f, x=%8.2f, y=%12.6f, z=%12.6f\n",
muzzleSpeed, s.t, s.x, s.y, s.z);
}
// *** Syntax description for findAngle() function
void findAngleUsage() {
fprintf(stderr, "findAngle usage:\n");
fprintf(stderr, "findAngle(degrees, \"degrees\", direction, \"of\", direction)\n");
fprintf(stderr, "possible combinations are:\n");
fprintf(stderr, "\"north\" \"of\" \"east\" \n");
fprintf(stderr, "\"east\" \"of\" \"north\" \n");
fprintf(stderr, "\"west\" \"of\" \"north\" \n");
fprintf(stderr, "\"north\" \"of\" \"west\" \n");
fprintf(stderr, "\"south\" \"of\" \"west\" \n");
fprintf(stderr, "\"west\" \"of\" \"south\" \n");
fprintf(stderr, "\"east\" \"of\" \"south\" \n");
fprintf(stderr, "\"south\" \"of\" \"east\" \n");
fprintf(stderr, "\n\n");
exit(1);
}
// *** Converts an English description of a direction into a polar coordinate angle
double findAngle(
double degrees,
const char * degreesT, // Must be the literal "degrees"
const char * direction1,
const char * ofT, // Must be the literal "of"
const char * direction2
) {
if(strcasecmp(degreesT, "degrees")) findAngleUsage();
if(strcasecmp(ofT, "of")) findAngleUsage();
double angle = 0;
if ((!strcasecmp(direction1, "north"))&&(!strcasecmp(direction2, "east"))) {angle = degrees;}
else if((!strcasecmp(direction1, "east"))&&(!strcasecmp(direction2, "north"))) {angle = 90-degrees;}
else if((!strcasecmp(direction1, "west"))&&(!strcasecmp(direction2, "north"))) {angle = 90+degrees;}
else if((!strcasecmp(direction1, "north"))&&(!strcasecmp(direction2, "west"))) {angle = 180-degrees;}
else if((!strcasecmp(direction1, "south"))&&(!strcasecmp(direction2, "west"))) {angle = 180+degrees;}
else if((!strcasecmp(direction1, "west"))&&(!strcasecmp(direction2, "south"))) {angle = 270-degrees;}
else if((!strcasecmp(direction1, "east"))&&(!strcasecmp(direction2, "south"))) {angle = 270+degrees;}
else if((!strcasecmp(direction1, "south"))&&(!strcasecmp(direction2, "east"))) {angle = 360-degrees;}
else {findAngleUsage();}
return(angle);
}
// *** Test jig for the findAngle() function
void testAngles(void) {
printf("%5.2f\n", findAngle(15, "degrees", "north", "of", "east"));
printf("%5.2f\n", findAngle(15, "degrees", "east", "of", "north"));
printf("%5.2f\n", findAngle(15, "degrees", "west", "of", "north"));
printf("%5.2f\n", findAngle(15, "degrees", "north", "of", "west"));
printf("%5.2f\n", findAngle(15, "degrees", "south", "of", "west"));
printf("%5.2f\n", findAngle(15, "degrees", "west", "of", "south"));
printf("%5.2f\n", findAngle(15, "degrees", "east", "of", "south"));
printf("%5.2f\n", findAngle(15, "degrees", "south", "of", "east"));
printf("%5.2f\n", findAngle(15, "degrees", "frunk", "of", "east"));
}
// *** This is the cannon aimer for lower trajectories.
// *** It works by repeatedly calculating the error, re-aiming accordingly, and taking another shot
const Vector planShot_low(
const Vector target, // vector location of the target
const double muzzleSpeed // scalar speed of the cannon shell
) {
// Declare and init vars
Vector hitPoint(0,0,0); // Where the shell crosses the target height going down
Vector prevHitPoint(0,0,0); // Hit point on previous iteration (break logic)
double xyAngleHitPoint; // xy angular direction of hitpoint
double xyAngleDeviation; // xy angular deviation between target and hitpoint
Vector muzzleVelocity(0,0,0); // Vector muzzle velocity producing the current hit point
double zAngle = -89; // elevation angle, start very low for shooting into valley
double xyAngle; // directional angle
double closeness = 100000000; // scalar distance between hit point and target,
// start huge to start loop
double overshoot = -1; // difference between shot distance and target distance,
// start with undershoot
double prevOvershoot; // previous iteration's overshoot
double interval = 8; // amount to increase elevation angle each iteration
double intervalFactor = 3.0; // divisor to decrease interval on new overshoots or undershoots
int maxIterations = 2000; // point to assume you're infinitely looping
int iteration; // iteration number
int blownBack; // flag -- was it blown by high wind backward away from target
// Find angular location of target
double xyAngleTarget = atan2(target.gety(), target.getx()) * Conversions::degreesPerRadian();
xyAngle = xyAngleTarget;
// Repeatedly estimate and shoot until you hit within 1 cm (.01) of target
iteration=0;
while(closeness > .01) {
// set previous values to current
prevHitPoint = Vector(hitPoint.getx(), hitPoint.gety(), hitPoint.getz());
prevOvershoot = overshoot;
// For this elevation angle, make sure the xy angular deviation is less/eq the interval/10
xyAngleDeviation = 99999;
while(xyAngleDeviation > fabs(interval)/10.0) {
if(iteration++ > maxIterations) break;
// printf("diag %d: zAngle=%10.6f, xyAngle=%10.6f, blBk=%d, closeness=%9.5f\n",
// iteration, zAngle, xyAngle, blownBack, closeness);
// Translate angles to a velocity vector
double zt = muzzleSpeed * sin(zAngle/Conversions::degreesPerRadian());
double xyt = muzzleSpeed * cos(zAngle/Conversions::degreesPerRadian());
double xt = xyt * cos(xyAngle/Conversions::degreesPerRadian());
double yt = xyt * sin(xyAngle/Conversions::degreesPerRadian());
muzzleVelocity = Vector(xt, yt, zt);
// Take the shot
struct Snapshot s = shootCannon(Vector(0,0,0), muzzleVelocity, target.getz(), 0.1, 0);
hitPoint = Vector(s.x, s.y, s.z);
// Calculate the xy angle of deviation between hitpoint and target
xyAngleDeviation = 0;
double xyAngleCorrection = 0;
xyAngleHitPoint = xyAngleTarget;
blownBack = 0;
// don't correct if shot never achieved altitude
if(hitPoint.getMagnitude() > 0.00000000001) {
xyAngleHitPoint = atan2(hitPoint.gety(), hitPoint.getx()) * Conversions::degreesPerRadian();
xyAngleDeviation = xyAngleHitPoint - xyAngleTarget;
xyAngleCorrection = xyAngleDeviation/2;
if((xyAngleDeviation > 90) || (xyAngleDeviation < -90)) blownBack = 1;
}
// Adjust xy angle for next iteration
if(!blownBack) xyAngle -= xyAngleCorrection;
}
//Below this point work with elevation angle incrementation and loop termination
// Determine how close you got (scalar)
closeness = target.scalarDistanceFrom(hitPoint);
// Calculate how much short or long
double tempAngle = (xyAngleHitPoint - xyAngleTarget)/Conversions::degreesPerRadian();
overshoot = hitPoint.getMagnitude() * cos(tempAngle) - target.getMagnitude();
// Reverse and shorten interval on over to undershoot transition or vice versa
if(prevOvershoot * overshoot < 0) {
interval /= intervalFactor;
interval *= -1;
}
// Increment zAngle
zAngle += interval;
// The rest of the loop detects and acts on loop termination criteria
// Detect total inability to achieve the target at any trajectory
if((zAngle > 180) || (zAngle < -180)) {
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
// Detect and brake out of infinite loop
if(iteration > maxIterations) {
fprintf(stderr, "Program error: Hung while trying to calculate shot for\n");
fprintf(stderr, "target (%8.3f, %8.3f, %8.3f).\n",
target.getx(), target.gety(), target.getz());
fprintf(stderr, "Call programmer: internal error\n\n");
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
}
return(muzzleVelocity);
}
// *** This is the cannon aimer for higher trajectories.
// *** It works by repeatedly calculating the error, re-aiming accordingly, and taking another shot
// *** Arg muzzleVelocityA must be the muzzle velocity that produces the low trajectory hit.
// *** Therefore this function can only be called after a call to planShot_low().
const Vector planShot_high(
const Vector target, // vector location of the target
const Vector muzzleVelocityA // VECTOR speed of the cannon shell
) {
// Guarantee overshoot
struct Snapshot s;
Vector muzzleVelocity = muzzleVelocityA;
double muzzleSpeed = muzzleVelocity.getMagnitude();
double desiredCloseness = .01; // stop when you're this close
double closeness = 0; // scalar distance between hit point and target
double prevCloseness;
while(closeness <= desiredCloseness) {
muzzleVelocity.addVector(Vector(0,0,1));
s = shootCannon(Vector(0,0,0), muzzleVelocity, target.getz(), 0.1, 0);
closeness = target.scalarDistanceFrom(Vector(s.x, s.y, s.z));
}
double overshoot = 1; // difference between shot distance and target distance
// Calculate some values
Vector hitPoint = Vector(s.x, s.y, s.z);
double xyAngle = atan2(muzzleVelocity.gety(), muzzleVelocity.getx()) * Conversions::degreesPerRadian();
double xySpeedTemp = sqrt(muzzleVelocity.getx() * muzzleVelocity.getx() +
muzzleVelocity.gety() * muzzleVelocity.gety());
double zAngle = atan2(muzzleVelocity.getz(), xySpeedTemp) * Conversions::degreesPerRadian();
double xyAngleTarget = atan2(target.gety(), target.getx()) * Conversions::degreesPerRadian();
// Declare and init other vars
Vector prevHitPoint(s.x, s.y, s.z); // Hit point on previous iteration (break logic)
double prevOvershoot; // previous iteration's overshoot
double interval = 1; // amount to increase elevation angle each iteration
double intervalFactor = 3.0; // divisor to decrease interval on new over or undershoots
int maxIterations = 2000; // point to assume you're infinitely looping
int iteration; // iteration number
int blownBack = 0; // Blown behind launch.
double xyAngleDeviation; // xy angular deviation between target and hitpoint
double xyAngleHitPoint; // angle of the hit point (for xy angle corrections)
// Repeatedly estimate and shoot until you hit within 1 cm (.01) of target
iteration=0;
while(closeness > desiredCloseness) {
xyAngleDeviation = 999999;
while(xyAngleDeviation > fabs(interval)/10) {
if(iteration++ > maxIterations) break;
// printf("diag %4d: zAngle=%10.6f, xyAngle=%10.6f, closeness=%9.5f, blBk=%d\n",
// iteration, zAngle, xyAngle, closeness, blownBack);
// set previous values to current
prevHitPoint = Vector(hitPoint.getx(), hitPoint.gety(), hitPoint.getz());
prevOvershoot = overshoot;
prevCloseness = closeness;
// Take the shot
struct Snapshot s = shootCannon(Vector(0,0,0), muzzleVelocity, target.getz(), 0.1, 0);
hitPoint = Vector(s.x, s.y, s.z);
// Calculate the xy angle of deviation between hitpoint and target
xyAngleHitPoint = xyAngleTarget;
xyAngleDeviation = 0;
double xyAngleCorrection = 0;
blownBack = 0;
// Correct left-right if shot achieved altitude
if(hitPoint.getMagnitude() > -0.00000000001) {
xyAngleHitPoint = atan2(hitPoint.gety(), hitPoint.getx()) * Conversions::degreesPerRadian();
xyAngleDeviation = xyAngleHitPoint - xyAngleTarget;
xyAngleCorrection = xyAngleDeviation/2;
if((xyAngleDeviation > 90) || (xyAngleDeviation < -90)) blownBack = 1;
}
// Adjust xy angles for next iteration
if(!blownBack) xyAngle -= xyAngleCorrection;
// Translate angles to a new velocity vector
double zt = muzzleSpeed * sin(zAngle/Conversions::degreesPerRadian());
double xyt = muzzleSpeed * cos(zAngle/Conversions::degreesPerRadian());
double xt = xyt * cos(xyAngle/Conversions::degreesPerRadian());
double yt = xyt * sin(xyAngle/Conversions::degreesPerRadian());
muzzleVelocity = Vector(xt, yt, zt);
// Determine how close you got (scalar)
closeness = target.scalarDistanceFrom(hitPoint);
}
// Calculate how much short or long
if(blownBack) {
overshoot = -1; // If it's blown back, it's definitely an undershoot
}
else
{
double tempAngle = (xyAngleHitPoint - xyAngleTarget)/Conversions::degreesPerRadian();
overshoot = hitPoint.getMagnitude() * cos(tempAngle) - target.getMagnitude();
}
if(prevOvershoot * overshoot < 0) { // if switched between over and under or vice versa
interval /= intervalFactor;
interval *= -1;
}
// The rest of the loop detects and acts on loop termination criteria
zAngle += interval; // raise elevation angle
// Detect total inability to achieve the target at any trajectory
if((zAngle > 180) || (zAngle < -180)) {
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
// Detect and brake out of infinite loop
if(iteration > maxIterations) {
fprintf(stderr, "Program error: Hung while trying to calculate shot for\n");
fprintf(stderr, "target (%8.3f, %8.3f, %8.3f).\n",
target.getx(), target.gety(), target.getz());
fprintf(stderr, "Call programmer: internal error\n\n");
closeness=0; //terminate loop
muzzleVelocity=Vector(0,0,0); //reflect inability to hit
}
}
return(muzzleVelocity);
}
// *** This is a test jig for the planShot() function
void attempt(
const double muzzleSpeed, //scalar speed of shell on firing
const double degrees,
const char * degreesT, // literal "degrees"
const char * direction1,
const char * ofT, // literal "of"
const char * direction2,
const double targetHeight // height of target
) {
double targetAngle = findAngle(degrees, "degrees", direction1, "of", direction2);
printf("\n\nangle is %f\n\n", targetAngle);
Vector target(
cos(targetAngle/Conversions::degreesPerRadian()),
sin(targetAngle/Conversions::degreesPerRadian()),
1
);
target.timesScalar(1.4142135 * Conversions::kilometersPerMile()* Conversions::metersPerKilometer());
target.setz(targetHeight);
Vector muzzleVelocity = planShot_low(target, muzzleSpeed);
struct Snapshot s= shootCannon(Vector(0,0,0), muzzleVelocity, targetHeight, 0.1, 0);
double windDirection = atan2(wi.getWindVelocity().gety(), wi.getWindVelocity().getx());
windDirection *= Conversions::degreesPerRadian();
printf("\nWind velocity=%9.4f at %f degrees\n", wi.getWindVelocity().getMagnitude(), windDirection);
printf(" Proof: Low shot... t=%8.3f x=%8.3f y=%8.3f z=%8.3f\n",s.t,s.x,s.y, s.z);
printf("Muzzle velocity=(%7.2f, %7.2f, %7.2f), zAngle=%8.3f\n",
muzzleVelocity.getx(),
muzzleVelocity.gety(),
muzzleVelocity.getz(),
Conversions::degreesPerRadian() * atan2(muzzleVelocity.getz(),
sqrt(muzzleVelocity.getx()*muzzleVelocity.getx() + muzzleVelocity.gety()*muzzleVelocity.gety())));
printf("\n-----------------------------------------------------------------\n");
muzzleVelocity = planShot_high(target, muzzleVelocity);
s= shootCannon(Vector(0,0,0), muzzleVelocity, targetHeight, 0.1, 0);
printf("Proof High shot... t=%8.3f x=%8.3f y=%8.3f z=%8.3f\n",s.t,s.x,s.y, s.z);
printf("Muzzle velocity=(%7.2f, %7.2f, %7.2f), zAngle=%8.3f\n",
muzzleVelocity.getx(),
muzzleVelocity.gety(),
muzzleVelocity.getz(),
Conversions::degreesPerRadian() * atan2(muzzleVelocity.getz(),
sqrt(muzzleVelocity.getx()*muzzleVelocity.getx() + muzzleVelocity.gety()*muzzleVelocity.gety())));
}
void main() {
wi.setWindVelocity(Vector(0,0,0));
double muzzleSpeed = getMuzzleSpeed();
confirmMuzzleSpeed(muzzleSpeed);
wi.setWindVelocity(Vector(14.14, -14.14, 0));
attempt(muzzleSpeed, 0, "degrees", "north", "of", "east", 777);
attempt(muzzleSpeed, 45, "degrees", "north", "of", "east", 777);
attempt(muzzleSpeed, 90, "degrees", "north", "of", "east", 777);
attempt(muzzleSpeed, 45, "degrees", "north", "of", "west", 777);
attempt(muzzleSpeed, 0, "degrees", "north", "of", "west", 777);
attempt(muzzleSpeed, 45, "degrees", "south", "of", "west", 777);
attempt(muzzleSpeed, 90, "degrees", "south", "of", "east", 777);
attempt(muzzleSpeed, 45, "degrees", "south", "of", "east", 777);
printf("\n\n");
}
|
| Language | Features |
| C | Low level, high performance language for systems and performance application programming. |
| C++ | A "better C than C", that sacrifices a little performance for object orientation and other features. |
| Perl | A highly productive interpreter with much better performance than you'd expect from an interpreter. Perl is installed by default on most Linux, Unix, and BSD computers. |
| Python | Another highly productive interpreter with high performance. Much more readable than Perl, useful for larger applications, but not as ubiquitous as Perl. |
| TCL/TK | Graphical output -- great for forms and the like. You can get the best of both worlds by using TCL/TK as the screen building component for Perl or Python. |
| wxPython | A graphic library for building GUI python applications. Very rapid development. |
| DBI::DBD | This sometimes doesn't come with the Linux distribution, but it's easily downloadable. It enables Perl to easily interface with most databases, including PostgreSQL and MySQL. |
| PostgreSQL | A very complete SQL DBMS implementation. |
| MySQL | A fast and lightweight SQL DBMS implementation. |
| PHP | An extremely productive language for creating web applications that interface to databases. It can also be used at the command line as a scripting language. |
| Java | Sun's Java -- the write once run anywhere language for web programming and other programming. |
| Ruby | A fully OOP language able to create apps similar to those created by Python or Perl. |
There are many specialized languages and tools for writing TK apps and other things.
The bottom line is this. If you get an idea and want to start experimenting it immediately, a Linux box is what you want. Your Linux box has many languages and libraries, all installed, for zero dollars.
-:-:-
My whole life I've bought records, and later CD's, on a regular basis. Once I got a VCR, I bought film videos on a regular basis. Whenever I went to Best Buy or Walmart, I always cruised the CD isles looking for a kewl 90's anthology, or house music, or 80's stuff, or whatever. At Walmart I always went down the video isle looking for a kewl new horror movie.
That all came to a halt several months ago when I boycotted everything from the music publishers, the film studios, and Disney.
These days I don't go near those isles, because I can't buy anything there. I'm boycotting the guys who are trying to destroy Linux. Six months ago it was called SSSCA, today it's called CBDTPA, but it's all the same -- it makes Linux and most Open Source illegal. It locks up Microsoft's monopoly tighter than it's ever been. It's not about piracy, it's about control. It's about the ultimate selfishness -- the willingness to destroy the works of innocent third parties just to maximize wealth.
Now, when my kids borrow my CD's, I remind them to take good care of them, because if my CD's are scratched I can't replace them. My wife still buys videos for the kids, but I don't. And I used to be their major source. So like it or not, the kids suffer. But when I told them the reason for my boycott, my 9 year old daughter said incredulously, "they're trying to outlaw Linux? That's terrible!". The kids understand my boycott.
-:-:-
As I mentioned, this atrocity has been foisted on American society by
the greedmongers in the record industry, the movie industry, and Disney.
But who were their congressional accomplices? The main guy is Fritz Hollings
from South Carolina. Here' s a list of the people killing Open Source for
the benefit of the greedy:
| Senator | - | Party | - | State |
| Fritz Hollings | D | South Carolina | ||
| Ted Stevens | R | Alaska | ||
| Daniel Inouye | D | Hawaii | ||
| John Breaux | D | Lousiana | ||
| Bill Nelson | D | Florida | ||
| Dianne Feinstein | D | California | ||
| Mr. HOLLINGS (for himself, Mr. STEVENS, Mr. INOUYE, Mr. BREAUX, Mr. NELSON of Florida, and Mrs. FEINSTEIN) introduced the following bill, which was read twice and referred to the Committee on |
If you live in one of these states, work to unseat these people. I'll certainly take a strong stand in campaigning against Nelson in his next run.
Hollings is the CBDTPA lead, and what a shame for South Carolina. South Carolina's governor, Jim Hodges, is involved in several lawsuits against the federal government because the feds are trying to transport nuclear waste from Colorado to South Carolina. According to one article (see URL's section), Mr. Hodges promises to physically prevent the feds from dumping this nuclear waste in his state.
How could a state have so little clout that it becomes the nation's nuclear dumpsite? Perhaps they have weak senators. Or perhaps their junior senator spent his time and political capital pandering to California fat cats instead of concentrating on his state's need to keep out radioactive waste. All of you from South Carolina, ponder that during the next election.
And when Hollings is replaced by a strong senator who cares more about his home state than he does about the L.A. fat cats, perhaps the feds will cast an eye on Florida as the nation's recipient of of nuclear waste. That's one reason I will be supporting Bill Nelson's opponent in the next election, whomever that opponent may be. I usually vote Democratic, but I'll work very hard for the election of Nelson's Republican opponent.
If you live in California, and work in technology, or just enjoy a robust economy, you're being sold down the river by your senator, Dianne Feinstein. CBDTPA is a disaster for the technology industry. By her support of CBDTPA, Feinstein is taking money out of your pocket and giving it to the movie studios and record companies. If you live in Silicon Valley and owe more on your house than its current market value, consider how Feinstein's actions will exacerbate your plight. Even if you live in Los Angeles, you're probably about to be shafted by Feinstein. Los Angeles has a huge technology sector that will be screwed by CBDTPA.
I see no end to my boycott. Even if CBDTPA is defeated, the greedmongers will try again. They're now trying to legislate copy protection into Digital to Analog converter chips (see URL's section). The boycott is difficult on me, and you may not wish to join me. But you can campaign, and you can vote. You can show these politicians that they can't shaft the people and keep their jobs. Prove that to them, and you'll finally have the representative government you deserve.
Hollings, Stevens, Inouye, Breaux, Nelson, and Feinstein. Do the right thing -- throw the bums out!
Any article submitted to Troubleshooting Professional Magazine must be licensed with the Open Publication License, which you can view at http://opencontent.org/openpub/. At your option you may elect the option to prohibit substantive modifications. However, in order to publish your article in Troubleshooting Professional Magazine, you must decline the option to prohibit commercial use, because Troubleshooting Professional Magazine is a commercial publication.
Obviously, you must be the copyright holder and must be legally able to so license the article. We do not currently pay for articles.
Troubleshooters.Com reserves the right to edit any submission for clarity or brevity, within the scope of the Open Publication License. If you elect to prohibit substantive modifications, we may elect to place editors notes outside of your material, or reject the submission, or send it back for modification. Any published article will include a two sentence description of the author, a hypertext link to his or her email, and a phone number if desired. Upon request, we will include a hypertext link, at the end of the magazine issue, to the author's website, providing that website meets the Troubleshooters.Com criteria for links and that the author's website first links to Troubleshooters.Com. Authors: please understand we can't place hyperlinks inside articles. If we did, only the first article would be read, and we can't place every article first.
Submissions should be emailed to Steve Litt's email address, with subject line Article Submission. The first paragraph of your message should read as follows (unless other arrangements are previously made in writing):
Copyright (c) 2001 by <your name>. This material may be distributed only subject to the terms and conditions set forth in the Open Publication License, version Draft v1.0, 8 June 1999 (Available at http://www.troubleshooters.com/openpub04.txt/ (wordwrapped for readability at http://www.troubleshooters.com/openpub04_wrapped.txt). The latest version is presently available at http://www.opencontent.org/openpub/).
Open Publication License Option A [ is | is not] elected, so this document [may | may not] be modified. Option B is not elected, so this material may be published for commercial purposes.
After that paragraph, write the title, text of the article, and a two sentence description of the author.