BIG BALLISTICS
Engine Background

Exterior Ballistics with MICROCOMPUTERS by W.J. Jurens

ROY L. NAVY HAS A THEORY concerning the loss of HMS Hood, but cannot confirm it without knowing the ballistics of the Bismarck's rifles. Jack Numbercruncher, having carefully studied Nathan Okun's now classic articles on armor penetration, is ready to compute the immunity zones for Austro-Hungarian battleships, but lacks ballistic information. Larry Longdrop wants to estimate the effects of aerial bombing on various ship designs, but is unfamiliar with the computational methods to employ. Despite the fact that gun systems lately have been almost entirely replaced by the missile as the primary weapons system of most surface warships, the fact remains that the naval rifle was considered to be one of the major arbiters of ship to ship combat for the majority of the time period of interest to readers of this journal. It is therefore somewhat surprising that until now there has been virtually nothing published concerning the detailed characteristics of these vital components of the naval arsenal. In spite of certain recent and commendable efforts to shed at least tabular light on the subject, it is still true that specific details on the armament systems installed aboard the ships being considered have been either heavily glossed or completely ignored in recent publications on ship design and construction.

The grand problem of the ballistician can be simply stated as follows: given a projectile of known physical properties such as size, shape, and mass, projected from a rifle at a known initial velocity (the determination of which is a problem in INTERIOR ballistics), the task is to compute the remaining velocity, angle of fall, time of flight, and range at any point on the resulting trajectory under specified atmospheric conditions. In the absence of other information, the official Range table produced by the country in question must be considered the definitive definition of the ballistics of any gun weapons system. Sadly, in a large percentage of cases, official range tables remain classified, inaccurate, or frankly non-existent. Although it has always been possible in principle to recompute or create from scratch those tables that were otherwise unavailable, it was only with the advent of the small, reasonably priced microcomputer that such a path became practically feasible. This paper is intended to present what is believed to be the first generally available microcomputer program for exterior ballistics investigations in naval history, accompanied by a sufficient theoretical base to allow the prospective user to generate accurate and complete range tables for any reasonable gun weapon system he or she might wish to study.

Source code for MASTER EXTERIOR BALLISTICS program in AppleSoft BASIC

(Download the source for the original plus a modified version with more drag tables)

10 HOME : VTAB 3: HTAB 4: INVERSE : PRINT "MASTER EXTERIOR BALLISTICS PROGRAM"
20 NORMAL
30 PRINT "*************************************"
40 PRINT "COPYRIGHT (C) (1983) W. J. JURENS"
50 PRINT " 62 FIDLER AVENUE"
60 PRINT " WINNIPEG, MANITOBA, CANADA"
70 PRINT " R3J 2R7"
80 PRINT " PH. 204-837-3125"
90 PRINT "*************************************"
100 PRINT
110 PRINT " PRESS ANY KEY TO LOAD THE PROGRAM": PRINT
120 GET A$
130 DIM M(50): DIM K(50): REM DIMENSION VARIABLES FOR M=MACH NO. AND K= DRAG COEFFICIENT FROM GRAPH OF AMERICAN KD8 FUNCTION
140 M(0) = 0.3:K(0) = 0.0825:M(1) = 0.5:K(1) = 0.0825:M(2) = 0.7:K(2) = 0.0825:M(3) = 0.8:K(3) = 0.0825:M(4) = 0.9:K(4) = 0.0825:M(5) = 0.95:K(5) = 0.094
150 M(6) = 1:K(6) = 0.140:M(7) = 1.05:K(7) = 0.172:M(8) = 1.1:K(8) = 0.175:M(9) = 1.15:K(9) = 0.175:M(10) = 1.2:K(10) = 0.1725:M(11) = 1.25:K(11) = 0.1690
160 M(12) = 1.30:K(12) = 0.166:M(13) = 1.35:K(13) = 0.1625:M(14) = 1.40:K(14) = 0.1605:M(15) = 1.45:K(15)=0.157
170 M(16) = 1.5:K(16) = 0.1535:M(17) = 1.55:K(17) = 0.151:M(18) = 1.6:K(18) = 0.149:M(19) = 1.65:K(19) = 0.146:M(20) = 1.7:K(20) = 0.143:M(21) = 1.75:K(21) = 0.14
180 M(22) = 1.8:K(22) = 0.138:M(23) = 1.85:K(23) = 0.136:M(24) = 1.9:K(24) = 0.1335:M(25) = 1.95:K(25) = 0.132:M(26) = 2:K(26) = 0.129:M(27) = 2.05:K(27) = 0.127:M(28) = 2.10:K(28) = 0.125:M(29) = 2.15:K(29) = 0.123:M(30) = 2.2:K(30) = 0.121
190 M(31) = 2.25:K(31) = 0.119:M(32) = 2.3:K(32) = 0.1175:M(33) = 2.35:K(33) = 0.116:M(34) = 2.4:K(34) = 0.1145:M(35) = 2.45:K(35)= 0.112:M(36) = 2.5:K(36) = 0.111
200 M(41) = 2.75:K(41) = 0.104:M(42) = 2.8:K(42) = 0.103:M(43) = 3.25:K(43) = 0.09:M(44) = 3.8:K(44) = 0.08:M(45) = 4.15:K(45) = 0.075:M(46) = 4.40:K(46) = 0.0735
210 PRINT
220 RD = 57.295779513: REM DEG/RAD CONVERSION
230 DF = 1.2250: REM DENSITY IN KG/M^3
240 Z4 = 1.34279408E - 18:Z3 = - 9.87941429E - 14:Z2 = 3.90848966E - 9:Z1 = - 9.69888125E - 5: REM POLYNOMIAL TERMS FOR ICAO STD ATMOSPHERE
250 INPUT " INPUT INITIAL VELOCITY (M/SEC) "; VO
260 PRINT
270 INPUT " INPUT ANGLE OF DEPARTURE (DEG) "; LD
280 PRINT
290 INPUT " INPUT PROJ. CALIBER (MM) "; DS
300 D = DS / 10: REM CAL MUST BE IN CENTIMETERS FOR FURTHER CALCS
310 PRINT
320 INPUT " INPUT PROJ. MASS (KG) ";M
330 PRINT
340 INPUT " INPUT FORM FACTOR ";FF
350 PRINT
360 INPUT " INPUT INTEGRATION INTERVAL (SEC) ";I
370 PRINT
380 INPUT " AUTO-SET INTERVALS? (Y/N) ";IC$: PRINT
390 IF IC$ = "Y" THEN INPUT " INPUT ERROR TOLERANCE (M/SEC) ";ET: PRINT
400 INPUT " INPUT AIR DENSITY FACTOR ";AD
410 PRINT
420 INPUT " INPUT ALTITUDE (METERS) ";AO
430 PRINT
440 PRINT " INPUT CHOICE OF DENSITY FUNCTIONS:"
450 PRINT " US. PRE-1945 STD= 1"
460 PRINT" BRITISH STANDARD = 2"
470 PRINT " I.C.A.O STANDARD = 3"
480 INPUT PW
490 HOME : VTAB 3: HTAB 4: INVERSE : PRINT "SUMMARY OF INITIAL CONDITIONS"
500 NORMAL : PRINT
510 PRINT " INITIAL VELOCITY = ";VO;" METERS/SEC."
520 PRINT " ANGLE OF DEPARTURE = ";LD;" DEGREES"
530 PRINT " PROJECTILE CALIBER = ";DS;" MM"
540 PRINT " PROJECTILE WEIGHT = ";M;" KG"
550 PRINT " FORM FACTOR = ";FF
560 PRINT " AIR DENSITY FACTOR = ";AD
570 PRINT " INTEGRATION INTERVAL:";I;" SEC"
580 PRINT " INITIAL ALTITUDE = ";AO;" METERS"
590 PRINT " ";X$
600 IF PW = 1 THEN PRINT " U.S. PRE-1945 DENSITY FUNCTION"
610 IF PW = 2 THEN PRINT " BRITISH STD. DENSITY FUNCTION"
620 IF PW = 3 THEN PRINT " I.C.A.O. DENSITY FUNCTION"
630 PRINT : PRINT : INVERSE : PRINT "PRESS 'G' TO GO"
640 NORMAL : PRINT "(ANY OTHER KEY WILL RECYCLE)": PRINT
650 GET PP$
660 IF PP$ < > "G" THEN HOME : GOTO 180
670 HOME : VTAB 3: HTAB 4: INVERSE : PRINT "TABULAR INTEGRATION OF TRAJECTORY"
680 PRINT "TIME RANGE HEIGHT ANGLE VELOCITY"
690 PRINT "SEC. METERS METERS DEC M/SECOND"
700 PRINT TT TAB(8) INT (10 * RG + 0.5) / 10 TAB(16) INT (10 * AO + 0.5) / 10 TAB(24) INT (1000 * (LD) + 0.5) / 1000 TAB(33) INT (1000 * VO + 0.5) / 1000
710 REM :ABOVE LINE TO PRINT INITIAL CONDITIONS
720 LD = LD / RD: REM CONVERT DEGREES TO RADIANS FOR CALCS
730 REM **BEGIN NUMERICAL INTEGRATION PROCEDURE**
740 C = M / (FF * AD * D ^ 2): REM BALLISTIC COEFFICIENT
750 IF TT = 0 THEN GOTO 830
760 PRINT TT TAB(8) INT (10 * RG + 0.5) / 10 TAB(16) INT (10 * AO + 0.5) / 10 TAB(24) INT (1000 * (L3 * RD) + 0.5) / 1000 TAB(33) INT (1000 * VO + 0.5) / 1000
770 REM :ABOVE LINE TO PRINT TYPICAL CONDITIONS
780 P3 = P2:P2 = P1:P1 = RG: REM RANGE ARRAY
790 Q3 = Q2:Q2 = Q1:Q1 = VO: REM VELOCITY ARRAY
800 S3 = S2:S2 = S1:S1 = L3: REM ANGLE ARRAY
810 T3 = T2:T2 = T1:T1 = TT: REM TIME ARRAY
820 U3 = U2:U2 = U1:U1 = AO: REM ALTITIUDE ARRAY
830 IF AO < 0 THEN GOSUB 1840
840 YE = YE + 1
850 IF YE = 5 THEN PRINT :YE = 0
860 REM :ABOVE LINES TO FORMAT PRINTOUT
870 Y = AO: REM SET UP FOR DENSITY
880 G = 9.80665 - (0.0000030665 * AO): REM ACCEL OF GRAVITY
890 CS = 344 - (0.004 * AO): REM SOUND SPEED (M/SEC)
900 XO = VO * COS (LD): REM HORZ COMPONENT
910 YO = VO * SIN (LD): REM VERT COMPONENT
920 REM :COMPUTE RETARDATION
930 KK = VO / CS: REM MACH NO.
940 FOR X = 1 TO 50
950 IF KK = < M(X) THEN GOTO 970
960 NEXT X
970 KD = (KK - M(X - 1)) / (M(X) - M(X-1)) * (K(X) - K(X-1)) + K(X-1): REM LINEAR INTERP. ON DRAG FUNCTION ARRAY
980 RO = KD * (DF / 10000) * VO ^ 2: REM RETARDATION IN M/SEC^2
990 Y = A1 + AO
1000 IF PW = 1 THEN GOSUB 1770
1010 IF PW = 2 THEN GOSUB 1790
1020 IF PW = 3 THEN GOSUB 1810
1030 DO = LL
1040 EO = RO / (C / DO): REM ACTUAL RETARDATION
1050 HO = EO * COS (LD): REM 1ST GUESS AT HORIZ COMPONENT OF RETARDATION
1060 JO = (EO * SIN (LD) + G): REM 1ST GUESS AT VERT COMPONENT OF RETARDATION
1070 X1 = XO - (HO * I): REM 1ST PRED. OF HORIZ VEL
1080 Y1 = YO - (JO * I): REM 1ST PRO). OF VERT VEL
1090 V1 = SQR (X1 ^ 2 + Y1 ^ 2)
1100 L1 = ATN (Y1 / X1): REM 1ST PRED. OF FINAL ANGLE
1110 REM ** BEGIN SECOND PREDICTION **
1120 M1 = (YO + Y1) / 2: REM ESTIMATE OF MEAN VERT VELOCITY
1130 A1 = (M1 * I): REM EST. ALTITUDE AT END OF INTERVAL
1140 Y = A1 + AO: REM SET UP TO COMPUTE DENSITY FOR ALTITUDE
1150 IF PW = 1 THEN GOSUB 1770
1160 IF PW = 2 THEN GOSUB 1790
1170 IF PW = 3 THEN GOSUB 1810
1180 D1 = LL
1190 KK = V1 / CS: REM MACH NO.
1200 FOR X = 1 TO 50
1210 IF KK = < M(X) THEN GOTO 1230
1220 NEXT X
1230 KD = (KK - M(X - 1)) / (M(X) - M(X - 1)) * (K(X) - K(X - 1)) + K(X - 1): REM LINEAR INTERP. ON DRAG FUNCTION ARRAY
1240 R1 = KD * (DF / 10000) * V1 ^ 2: REM RETARDATION IN M/SEC^2
1250 E1 = R1 / (C / D1): REM ACTUAL RETARDATION
1260 H1 = E1 * COS (L1): REM 2ND GUESS AT HORIZ RET. COMPONENT
1270 J1 = (E1 * SIN (L1) + G): REM 2ND GUESS AT VERT RET. COMPONENT
1280 H2 = (HO + H1) / 2: REM 3RD EST OF HORIZ RET COMPONENT
1290 J2 = (JO + J1) / 2: REM 3RD EST OF VERT RET COMPONENT
1300 X2 = XO - (H2 * I): REM 2ND PRED OF HORIZ VEL
1310 Y2 = YO - (J2 * I): REM 2ND PRED OF VERT VEL
1320 V2 = SQR (X2 ^ 2 + Y2 ^ 2): REM 2ND PRED OF FINAL VELOCITY
1330 L2 = ATN (Y2 / X2): REM 2ND PRED OF FINAL ANGLE
1340 REM ** BEGIN THIRD PREDICTION **
1350 M2 = (YO + Y2) / 2: REM 2ND EST OF MEAN VERT VEL
1360 A2 = M2 * I: REM EST ALT AT END OF INTERVAL
1370 Y = A2 + AO: REM SET UP FOR DENSITY CALC
1380 IF PW = 1 THEN GOSUB 1770
1390 IF PW = 2 THEN GOSUB 1790
1400 IF PW = 3 THEN GOSUB 1810
1410 D2 = LL
1420 KK = V2 / CS: REM MACH NO
1430 FOR X = 1 TO 50
1440 IF KK = < M(X) THEN GOTO 1460
1450 NEXT X
1460 KD = (KK - M(X - 1)) / (M(X) - M(X - 1)) * (K(X) - K(X - 1)) + K(X - 1): REM LINEAR INTERP. ON DRAG FUNCTION ARRAY
1470 R2 = KD * (DF / 10000) * V2 ^ 2: REM RETARDATION IN M/SEC^2
1480 E2 = R2 / (C / D2): REM ACTUAL RETARDATION
1490 H3 = E2 * COS (L2): REM 4TH GUESS AT HORIZ RET COMPONENT
1500 J3 = (E2 * SIN (L2) + G): REM 4TH GUESS AT VERT RET. COMPONENT
1510 H4 = (HO + H3) / 2: REM 5TH EST OF HORIZ RET COMPONENT
1520 J4 = (JO + J3) / 2: REM 5TH EST OF VERT RET. COMPONENT
1530 X3 = XO - (H4 * I): REM 3RD PRED OF HORIZ VEL
1540 Y3 = YO - (J4 * I): REM 3RD PRED OF VERT VEL
1550 V3 = SQR (X3 ^ 2 + Y3 ^ 2): REM 3RD PRED OF FINAL VELOCITY
1560 L3 = ATN (Y3 / X3): REM 3RD PRED OF FINAL ANGLE
1570 M3 = (YO + Y3) / 2: REM 3RD EST OF MEAN VERT VEL
1580 M4 = (XO + X3) / 2: REM 3RD EST OF MEAN HORIZ VEL
1590 FH = M4 * I: REM FINAL HORIZ POSN
1600 FV = M3 * I: REM FINAL VERT POSN
1610 LI = I: REM PUT INTERVAL IN SCRATCH VARIABLE *BEFORE* CHANGING VIA AUTO SET
1620 REM :** END OF 3RD PREDICTION **
1630 REM :DO AUTO SET OF INTERVAL
1640 IF IC$ < > "Y" THEN GOTO 1680
1650 IF ABS (V3 - V2) < 0.5 * ET THEN I = I * 1.2
1660 IF ABS (V3 - V2) > 0.8 * ET THEN I = I * 0.8
1670 IF ABS (V3 - V2) > ET THEN INVERSE : PRINT "ERROR TOLERANCE EXCEEDED": NORMAL : I = I / 2: REM IF ERROR EXCEEDED WARN AND REDUCE INTERVAL
1680 REM :** SETUP TO LOOP BACK **
1690 IF VE > ME THEN ME = VE: REM LINE TO FIND LARGEST VEL ERROR
1700 VE = ABS (V2 - V3): REM ERROR IN FINAL VELOCITY PREDICTIONS
1710 VO = V3
1720 LD = L3
1730 AO = AO + FV: REM FINAL ALTITUDE
1740 RG = RG + FH: REM SUM TOTAL RANGE
1750 TT = INT (1000 * (TT + LI) + 0.5) / 1000: REM SUM INTEGRATION INTERVALS FOR TOTAL TIME & CORRECT INTERPRETER ROUNDOFF ERROR
1760 GOTO 740
1770 LL = 10 ^ - (0.000045 * Y) : REM DENSITY OF ATMOSPHERE, MODEL 1 US PRE 1945 STD
1780 RETURN
1790 LL = 0.1 ^ (0.141 * (Y / 3048)) : REM DENSITY OF ATMOSPHERE, MODEL 2 BRITISH STD
1800 RETURN
1810 LL = Z4 * Y ^ 4 + Z3 * Y ^ 3 + Z2 * Y ^ 2 + Z1 * Y + 1 : REM DENSITY OF ATMOSPHERE, ICAO STD
1820 RETURN
1830 REM BEGIN FINAL INTERPOLATION ROUTINE
1840 Y(1) = P3:Y(2) = P2:Y(3) = P1
1850 X(1) = U3:X(2) = U2:X(3) = U1
1860 GOSUB 2050
1870 RL = B
1880 Y(1) = Q3:Y(2) = Q2:Y(3) = Q1
1890 GOSUB 2050
1900 SL = B
1910 Y(1) = S3:Y(2) = S2:Y(3) = S1
1920 GOSUB 2050
1930 AL = B
1940 Y(1) = T3:Y(2) = T2:Y(3) = T1
1950 GOSUB 2050
1960 TL = B
1970 PRINT : PRINT : HTAB 4: INVERSE : PRINT "INTERPOLATED TERMINAL CONDITIONS": NORMAL : PRINT
1980 PRINT "RANGE = "; INT (RL);" METERS"
1990 PRINT "HEIGHT = 0 METERS"
2000 PRINT "VELOCITY = "; INT (10 * SL + 0.5) / 10;" METERS/SEC."
2010 PRINT "ANGLE OF FALL = "; AL * RD;" DEGREES"
2020 PRINT "TIME = "; INT (TL * 100 + 0.5) / 100;" SECONDS"
2030 PRINT : HTAB 10: PRINT "FINISHED WITH PROBLEM": END
2040 REM BEGIN 3 POINT ROUTINE
2050 P = 3
2060 B = 0
2070 A = 0
2080 FOR J = 1 TO P
2090 T = 1
2100 FOR I = 1 TO P
2110 IF I = J THEN 2130
2120 T = T * (A - X(I)) / (X(J) - X(I))
2130 NEXT I
2140 B = B + T * Y(J)
2150 NEXT J
2160 RETURN

KD8 Type Projectile Dimensions

AppleSoft BASIC Emulator (http://www.calormen.com/jsbasic/)

Very cool. Thanks.

I'm glad to report that the program, albeit in somewhat expanded and enhanced form, e.g. translated into different computer languages, is still in use today. It remains, so far as I know, the only program in the public domain, available at reasonable cost, which actually does long-range artillery ballistics accurately. (As you have probably noticed, although there are a LOT of other ballistics programs 'out there', most or all are concerned with small-arms ballistic issues.

I hope you find it useful and even interesting to use. There are two errors in the text. My middle initial is "J" not "R", and one of the drag functions -- I can't remember which one -- is actually reproduced twice at the back of the text. This is not a problem, as one can quite easily enter new drag functions as required using the basic format of the tables already provided.

Again, thanks for sending this along...

Bill Jurens