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:
////////////////////////////////////////////////// |
