THE MATHEMATICS OF FALLING BOMBS
Harold Brochmann

As any high school student of Physics knows, plotting the trajectory followed by a falling bomb is easily done with the formulas:

x = vt, and
y = yo - (gt^2)/2

We assume that horizontal speed remains constant and is not affected one way or another by the bomb's simultaneous vertical motion. The distance that the bomb travels horizontally is found from the product of the horizontal velocity and the elapsed time. The vertical velocity of objects in free fall increases constantly and is not affected by the object's horizontal motion. The altitude at any one time is found by subtracting the product of the acceleration due to gravity and the square of the time from the initial altitude. Using 9.8 metres per second per second for g and some arbitrary value for v it is an easy matter to calculate and plot coordinates of representative points along a hypothetical bomb's trajectory. The path is parabolic. Another demonstration of applied mathematics in action.

There is only one problem with this - its rubbish. The formulas do not take into account the reality of air resistance, and therefore the results are pure fiction. In the real world falling objects do not keep accelerating; they reach what is called terminal velocity, and then don't fall any faster. Anyone who has tried skydiving knows this. Similarly, the horizontal velocity also slows. The topic came up during a recently during a discussion with a friend of mine whose hobby is target shooting. The question of what the early computers were used for was raised. He knew something that I did not; that ENIAC was used mostly to calculate the trajectories of artillery shells and bombs - a fascinating field called ballistics.

I had always assumed that air resistance is a function of the square of the velocity It seemed reasonable at the time!. This turns out not to be so at all. I quote from Mechanics and Properties of Matter by Reginald Stephenson:

When an object moves through a medium such as air, it is subjected to a resistance force. This resistance force is a function of the velocity and the shape of the moving object. There is no simple relationship between the velocity and the force, but at very low velocities it appears that the force varies approximately as the velocity, whereas at higher velocities it appears to be approximately proportional to the square or even higher powers of the velocity. In any practical situation the relationship between the resistance force and the velocity has to be determined experimentally.

It turns out that my friend has made quite a study of ballistics. He had a list of the Gavre Drag Values. The Gavre Drag Values are a set of numbers representing the effects of air resistance on a standard projectile at different velocities. The standard projectile was defined as a 2 inch long, 1 inch diameter cylinder with one end pointed. The values are empirical and were obtained by working backwards from experimental results. If you plot them you get a graph that looks like this:

The vertical axis represents air resistance. The horizontal axis is speed. The graph appears to be close to linear at its upper and lower extremes. It is steepest in the range around the velocity of sound. Interesting.

Ballistics calculations that determine the path followed by a falling object in air are of a type called 'iterative'. You consider the cumulative effects of what happens during small slices of time. The velocity at the end of any small slice of time is found by taking the velocity at the beginning and multiplying by a coefficient which incorporates the effects of air drag during that interval. This new value becomes the input for the next time slice and so on. Very small time slices are needed for valid results.

My friend had taken the original equations and re-written them into a form that would be amenable to FORTRAN programming. Here they are:

To illustrate how the equations work, let us consider what happens during the first one second in the case of a bomb in the shape of a standard projectile dropped at an altitude of 500 m from an airplane travelling 100 m/s. The Gavre values for 0 and 100m/s are 0.014 and 0.03241 respectively. For a complete listing of Gavre values consult the Pascal program listing at the end of the article. The equations don't actually work precisely for time slices as long as one second.

NXvel = IXvel * e-Kx*dt 
      = 100 * e-0.03241*1
      = 96.81

If there was no air resistance, the horizontal speed would remain at 100 m/s. With air resistance the horizontal speed is less - 96.81m/s....

If there was no air resistance, the horizontal distance travelled during that first second would be 100m. With air resistance, only 95.26m are covered horizontally.

Without air resistance the downward velocity at the end of one second would be 9.8m/s. In this case the speed is less than that.

Without air resistance, the distance covered during the first second would be the average speed * time: (0 + 9.8)/2 *1 = 4.9m. This would place the bomb at an altitude of 500 - 4.9 = 495.1m. In this case the altitude is slightly greater.

I wrote a Pascal Program which determines the coordinates that describe the path of a bomb dropped from an airplane in level flight. The ballistic coefficient, the plane's speed and the altitude as well as the duration of the time slice can all be specified. The default values are 1, 100, 500, 0.01 respectively. The program prints out two tables, the first,without air resistance, looks like this:

TIME INCREMENT: 0.010 SECONDS 
BALLISTIC CONSTANT: 1.00  
Time  Xacc  Yacc  Xspeed  Yspeed   Xpos    Alt.  
 0.00  0.00  9.80  100.00   0.00    0.00  500.00  
 1.00  0.00  9.80  100.00   9.80  100.00  495.10  
 2.00  0.00  9.80  100.00  19.60  200.00  480.40  
 3.00  0.00  9.80  100.00  29.40  300.00  455.90  
 4.00  0.00  9.80  100.00  39.20  400.00  421.60  
 5.00  0.00  9.80  100.00  49.00  500.00  377.50  
 6.00  0.00  9.80  100.00  58.80  600.00  323.60  
 7.00  0.00  9.80  100.00  68.60  700.00  259.90  
 8.00  0.00  9.80  100.00  78.40  800.00  186.40  
 9.00  0.00  9.80  100.00  88.20  900.00  103.10 
10.00  0.00  9.80  100.00  98.00  1000.00  10.00
Last two positions:
   1010.00      0.15
   1011.00     -0.84
Intrapolated horizontal distance to impact point:
    1010.2

Although the time slice being used is 0.01 sec, only whole second results are being recorded to save space. The horizontal acceleration is 0 (zero) .The vertical acceleration is constant at 9.8 m/sec/sec . The horizontal velocity is constant. The vertical velocity increases by 9.8 m/sec for each 1 sec time interval. The last calculated position has a negative altitude. The impact point is determined from the last two positions.

Here are the results of a run with air drag taken into account:

TIME INCREMENT: 0.010 SECONDS
BALLISTIC CONSTANT: 1.000
 Time   Xacc   Yacc Xspeed Yspeed    Xpos    Alt.
 0.00  -3.24   9.80 100.00   0.00    0.00  500.00
 1.00  -2.53   9.65  97.43   9.73   98.68  495.12
 2.00  -2.45   9.49  94.94  19.30  194.83  480.60
 3.00  -2.37   9.30  92.53  28.69  288.53  456.59
 4.00  -2.30   9.10  90.20  37.90  379.87  423.29
 5.00  -2.23   8.88  87.93  46.89  468.90  380.88
 6.00  -2.16   8.50  85.74  55.56  555.71  329.63
 7.00  -2.10   8.27  83.61  63.95  640.36  269.87
 8.00  -2.03   8.04  81.54  72.11  722.91  201.84
 9.00  -1.97   7.81  79.54  80.04  803.42  125.77
10.00  -1.92   7.58  77.59  87.73  881.96   41.89
Last two positions:
    917.45      0.74
    918.21     -0.17
Intrapolated horizontal distance to impact point:
     918.1

We note that the horizontal and vertical accelerations decrease.

Here is what happened when we graph the last two columns for each run.

We see that with no air resistance the bomb hits the ground just past the 1000 m mark, and with air drag its not quite 100 m short of that. The program calculates these to be at 1010.2 m and 918.1 m respectively.

We now turn to an examination of the Pascal program. If you want to try this program yourself I suggest you download the complete source code listing by using the link at the end of the article.

Here is the pseudo program:

program bomb;
{var}
{procedure initialize}
{function Gavre}
{procedure speeds_and_positions_air}
{procedure speeds_and_positions_vacuum}
{procedure update_variables}
{function finished}
{procedure find_impact}
{*Main*}
{initialize}
{speeds_and_positions_air}
{speeds_and_positions_vacuum}
{end.}

We start with the variable declarations, followed by the initializing procedure which gives initial values to the variables and also prints the headings for columns where the results of the calculations will be written..

var
	ixp, 	{initial x position}
	iyp, 	{initial y position}
	ixs, 	{initial x speed}
	iys, 	{initial y speed}
	dt, 	{time increment}
	nxp, 	{new x position}
	nyp, 	{new y position}
	nxs, 	{new x speed}
	nys, 	{new y speed}
	mxs, 	{mean x speed}
	mys, 	{mean y speed}
	xp, yp, 	{last values}
	bc, 	{ballistic coefficient}
	time:	{elapsed}
	extended;
	Ga: array[0..17] of real;
 const
	g = 9.8;	{metres / sec / sec}
 procedure initialize;
 begin
	dt := 0.01;	{seconds}
	time := 0;	{seconds}
	ixp := 0;	{position 0 horizontally}
	iyp := 500;	{altitude in metres}
	ixs := 100;	{metres / second}
	iys := 0;	{level flight}
	bc := 1;	{coefficient}
	Ga[0] := 0.01400;
	Ga[1] := 0.01991;
	Ga[2] := 0.02609;
	Ga[3] := 0.03557;
	Ga[4] := 0.05910;
	Ga[5] := 0.10968;
	Ga[6] := 0.16191;
	Ga[7] := 0.19482;
	Ga[8] := 0.21832;
	Ga[9] := 0.23611;
	Ga[10] := 0.25053;
	Ga[11] := 0.26301;
	Ga[12] := 0.27415;
	Ga[13] := 0.28502;
	Ga[14] := 0.29568;
	Ga[15] := 0.30694;
	Ga[16] := 0.31865;
	Ga[17] := 0.33496;
	writeln('TIME INCREMENT: ', dt : 5 : 3, ' SECONDS');
	writeln('BALLISTIC CONSTANT: ', bc : 5 : 3);
	write('Time' : 5, 'Xacc' : 7, 'Yacc' : 7);
	write('Xspeed' : 7, 'Yspeed' : 7);
	writeln('Xpos' : 8, 'Alt.' : 8);
 end;

I think that all of these are reasonably self-explanatory with the exception of the array Ga. Ga[0] contains the Gavre value for 0 velocity. Ga[1] for 50 f(m,s), Ga[2] for 100 f(m,s) and so on. I interpolated these metric values from the original fps table that my friend had. We now turn to the program itself:

{* Main *}
begin
 initialize;
 while not finished do
	begin
	 speeds_and_positions_air;
	{speeds_and_positions_vacuum;}
	 update_variables;
	end;
 find_impact
end.

In English it says:

Note that there are two procedures for doing calculations, one for vacuum and one for air. One or the other is activated by placing squiggly braces around the unwanted procedure call.

 procedure speeds_and_positions_vacuum;
 begin
	nxs := ixs;
	nxp := ixp + (nxs * dt);
	nys := iys + (g * dt);
	mys := (iys + nys) / 2;
	nyp := iyp - (mys * dt);
 end;

In English it says:

procedure update_variables;
	var
	 xacc, yacc: real;
 begin
	xacc := (nxs - ixs) / dt;
	yacc := (nys - iys) / dt;
	if abs(time - round(time)) < 0.001 then
	 begin
		write(time : 5 : 2);
		write(xacc : 7 : 2, yacc : 7 : 2);
		write(ixs : 7 : 2, iys : 7 : 2);
		writeln(ixp : 8 : 2, iyp : 8 : 2);
	 end;
	time := time + dt;
	xp := ixp;
	yp := iyp;
	ixp := nxp;
	iyp := nyp;
	ixs := nxs;
	iys := nys;
 end;

The horizontal and vertical accelerations are calculated. These are not really needed, but why not include them for interest? The results of the calculations are written in the appropriate columns. Finally, the initial values for the next time slice are equated to the calculated values in preparation for the next time slice.

There is an if... then condition around the part that writes things to the screen. This is so that only the values corresponding to whole second times periods will be written. This condition can be removed and you will get many more data points.This slows things down quite bit.

 function finished: boolean;
 begin
	finished := false;
	if nyp < 0 then
	 finished := true;
 end;

This is a function that simply returns true or false depending on whether or not the altitude (nyp) has reached a negative value. The procedure which finds the impact point intrapolates values between the last two positions. We assume that the curve is close enough to being linear for this to be ok.

 procedure find_impact;
 begin
	writeln('Last two positions:');
	writeln(xp : 10 : 2, yp : 10 : 2);
	writeln(nxp : 10 : 2, nyp : 10 : 2);
	xp := (yp / (yp - nyp) * (nxp - xp)) + xp;
	writeln('Extrapolated horizontal distance to impact point:');
	writeln(xp : 10 : 1);
 end;

Now for the calculations when air drag is taken into account:

procedure speeds_and_positions_air;
	var
	 Kx, Ky,: extended;
 begin
	Kx := Gavre(ixs) / bc;
	Ky := Gavre(iys) / bc;
	nxs := ixs * exp(-1 * Kx * dt);
	nxp := ixp + (nxs / Kx) * (1 - exp(-1 * Kx * dt));
	nys := ((iys - (g / Ky)) * exp(-1 * Ky * dt)) + (g / Ky);
	nyp := iyp - ((1 / (Ky * Ky)) * (g + (Ky * nys)) * (1 - exp(-1 * Ky * dt)) - ((g * dt) / Ky));
 end;

Here we have declared 2 more variables. These are local to the procedure. Kx and Ky are obtained by dividing the corresponding Gavre values by the ballistic coefficient. The Gavre values are obtained by a function which I'll discuss next. When the function is called, the appropriate velocity is passed to it; so for example, Gavre(325) returns the Gavre value for a velocity of 350 f(m,s). The rest of this procedure simply applies the formulas discussed earlier. Actually, they aren't simple at all!

 function Gavre (speed: extended): extended;
	var
	 i, 	{utility counter}
	 mini, maxi	{bounds of i}
	 : integer;
	 ming, maxg	{bounds of Ga value}
	 : real;
 begin
	for i := 0 to 17 do
	 begin
		if (i * 50 > speed) then
		 leave;
	 end;
	mini := i - 1;
	maxi := i;
	ming := Ga[mini];
	maxg := Ga[maxi];
	Gavre := (speed / (maxi * 50)) * (Ga[maxi] - Ga[mini]) + Ga[mini];
 end;

To give a specific example: Suppose you wanted to find the Gavre value for 320 m/sec. Ga[6] contains the value for 300 m and Ga[7] contains the value for 350 m. Assuming there is sufficient linearity here, we can find the value for 320 m by:

Gavre = (320 / (7 * 50)) * (Ga[7] - Ga[6]) + Ga[6]
      = (320/350) * (0.19482 - 0.16191) + 0.16191
      = 0.191999...

Anyway, after some local variables are declared, a loop finds the upper index value (7 in this example), and assigns it to variable maxi, with mini being 1 less (6). The corresponding Gavre values are assigned to maxg and ming, and then the are not is performed.

Well, there it is. There are many things here for the interested to explore. After the program has been entered and is running as advertised one might try plotting comparative trajectories for ÔbombsÕ of lower ballistic constants. It would not be very difficult to alter the program so that one could determine terminal velocities for free falling objects. Another challenge would be to determine the angle at which a gun should be elevated in order to achieve maximum range. One would think that precise answers to such questions would be influenced by the density of the air at different altitudes, its humidity and temperature, and possibly the spin imparted to the projectile by the rifles in a gun barrel. Apparently, though, these factors have relatively small significance and the equations we have discussed produce reasonable results.

My friend wrote the original version of this program in Fortran way back when because his competition rifle had been sighted in for 100 yards, and he wanted to know by how much to adjust the sight for 500 metres! He is now retired. And yes; he is a very good shot.

For the complete Pascal source code listing, click here