BASIC programs for battleship fire control and antiaircraft gun fire control (also a space battleship
fire control program)
By Ian W. Douglas, July 2011. Appendix (see end of article) added January 2012 and revised March 2012.
Further revised December 2011: these further revisions involve the programs TICTAC1.BAS, BATCON1.BAS,
ATICTAC.BAS, SPACBAT.BAS, and AIRBAT.BAS . I rewrote sections of these programs to enable the programs
to better handle situations where the target ship or aircraft are following perfect 0 degree, 90
degree, 180 degree and 270 degree course azimuths, as measured by the (own) ship or antiaircraft gun
tracking instrument and computer.
Further revised March 2012: standardized the flow and logic between IOELLIP.BAS and IOELOID.BAS to make
the flow and logic of these two programs as similar as possible, wherever possible. Fixed some bugs
in these two programs (division by zero was sometimes possible in the earlier versions of these programs,
in cases where a target or point was on an axis or had zero as one of the coordinates). Added programs
GUNEOD1.BAS and EIDGUN2.BAS, together with explanatory notes for these programs; these programs determine
whether a stationary airborne target is within range of an antiaircraft gun, assuming that the only forces
acting on the shell are the gun's initial impulse and gravity.
Further revised July 2012: clarified function (via remarks in program listing) of some of the parts of
certain programs; parts of some of the programs convert a vector and starting point to an end point.
I did not explain this clearly in earlier versions of programs.
Further revised December 2012: added more discussion of the problem of whether the target is within
range of an antiaircraft gun; added program UPARTR.BAS at end of Appendix; made some comments in program
listings more consistent internally; and alerted the reader to weaknesses in UPART.BAS, AIRBAT.BAS, and
UPARTR.BAS.
I have written some programs in BASIC which deal with the problem of tracking an enemy ship from a
moving ship, and for battleship gun fire control. The listings for the programs are included below
on this page and file. I also wrote some programs (included on this page as well) which solve the
same problem, except with aircraft tracking and antiaircraft gun fire control. Some of the REM (remark
or comment) lines in the program listings refer to flowcharts I made in order to write the programs, and
I would be willing to send text descriptions of the flowcharts to those who are interested (my flowcharts
are hand drawn and hand written and take up several pages, so they would be difficult to send graphically).
Readers who are familiar with flowcharting and programming can reconstruct my flowcharts from my program
listings, and those who wish to can send me their guesses to me as to what my flowcharts look like. I
should be able to send back my comments as to what is right or wrong about your guesses about my flowcharts.
I have sprinkled many comments throughout the program listings so as to make understanding of the programs'
functions and logic relatively easy to understand, assuming that the reader has some knowledge of BASIC
programming. I have test run the programs successfully using the Microsoft GW-BASIC interpreter.
The programs use a cross between the Dreyer plot (relative plot) and the Pollen plot (true plot).
I used a modified true plot which assumes that our own ship's heading is always 0 degrees/360 degrees;
however, given this, our own ship moves against a fixed background in a plot of the battle. If our own
ship does not turn, then this is a good "Occam's razor" approach to the problem in which we assume a 2
dimensional Cartesian coordinate system for the two ships. The origin of the coordinate system I used
puts our own ship at (50 nautical miles, 50 nautical miles) at the start of the problem so as to avoid
using negative numbers to represent points in space occupied by the ships and projectile. The projectile
does move in three dimensions; however, since gun and target are in 2 dimensions, the problem of the
projectile's motion is--for intents and purposes inherent in ship-to-ship firing--more or less two
dimensional. My program uses, and takes into account, the arc of trajectory.
There are two programs I wrote which address the battleship fire control problem. The first I called
"TICTAC1.BAS" because of the acronym TICTAC I used, which stands for "time interval compensating tachymetry."
This program uses data from a range and bearing finder (assuming a flat sea, a ship which does not pitch,
yaw, or roll, etc.) and our own ship's speed (negative speeds accepted by the program for our own ship going
backwards) to determine the enemy ship's starting coordinates, course and speed, and end coordinates.
The program allows for our own ship's speed when calculating the enemy ship's speed, direction, and
absolute coordinate data.
The second program I called "BATCON1.BAS" because it is the battleship fire control program. It uses
data from TICTAC.BAS, together with own ship's speed (the program accepts negative speeds for a ship going
in reverse, with forward heading still assumed to be 360 degrees) and muzzle velocity of the gun, to obtain
a fire and hit solution. The program assumes a flat sea, no air resistance, no wind speed and direction,
no Coriolis effect, acceleration due to gravity of 32 feet per second per second, and no rotating earth.
If the target is forward of the beam of our own ship, and if our own ship is moving forwards, then the
program uses the following method. The program uses trial and error to find a point uprange from the
target ship to shoot at--uprange by a distance equal to the distance travelled by our own ship during the
shell's time of flight, and along a line parallel to our own ship's course. At the outset of the problem
we have a "chicken and egg" situation. How far uprange (on the sea's surface) from the target to aim for
(I called this point "point C" in the program) depends on the shell's time of flight, but the shell's time
of flight depends on where point C is. The program reiterates itself until there is convergence between
time of shell's flight and a time to traverse from point C to the enemy at our own ship's speed.
If the target is on, or abaft the beam of, our own ship, and our own ship is moving forwards, then point C
is downrange of the target. Generally speaking, though, the above method still applies.
To put it another way: let's say that one wishes to do a drive by shooting from a boat at a buoy. The boat
is travelling due north at 10 feet per second, and when the buoy is abeam of the boat's gun, it will be due
west of the boat and gun. At that point, the buoy will be 6840 feet off the port beam. The gun's muzzle
velocity is 800 feet per second, and to make the projectile go those 6840 feet over the water's surface, the
gun must be elevated to 20 degrees, and the projectile will take 17.1 seconds to reach the target once fired.
(I did all these calculations using formulae used in the subsequent programs on this web page.) This means
that one must fire the gun directly off the port beam (i.e. due west of the boat) when the boat is 17.1 seconds
away from reaching a point where the buoy is off the port beam. Point "C" is the point on the water's surface
that you fire at, and it is 10 feet per second * 17.1 seconds = 17.1 feet due south of the buoy. (Distance
= rate * time). Point "C" is also 6840 feet due west of the boat at the instant the gun is fired.
If necessary, I could scan a diagram I made of the problem and sent it to those who are interested as a
GIF file. Also I would encourage readers to try the programs and see if they give answers which make sense.
As I mentioned above, I also wrote programs to do aircraft tracking (ATICTAC.BAS) and antiaircraft gun fire
control (AIRBAT.BAS). I also wrote what I call a "space battleship fire control program" (SPACBAT.BAS,
listing at the end of this file). This is similar to the antiaircraft gun fire control program I wrote
(I should note here that it even takes as input data the results of my aircraft tracking program) and
generates an antiaircraft or anti spacecraft gun firing solution, assuming you are in deep space with no gravity.
There is a strange phenomenon I have observed when running BATCON1.BAS. Let's say that our own ship is heading
due north, and that the target ship is behind us, abaft the port beam, and heading due west. The target ship
has a red flag on its bow and a blue flag on its stern. We have crossed the target ship's path and it is re-
ceding into the distance abaft the port beam. We look at the centreline of our own wake, and this is a line
extending due south into the distance. We sweep our gaze to the right (westward) and see the target ship.
Our own ship is travelling at 20 knots and the target ship is travelling at 5 knots. Due to the changing
relative positions of the target ship and own ship, and the geometry of the overall layout of the situation,
we see the target ship appear to get smaller; but it also appears to go backwards, apparently moving towards
the centreline of our own wake (you can verify this the next time you are out boating). Running BATCON1.BAS,
though, I noticed that in order to get a firing solution, one must make a further step. You have to aim at
a point apparently behind the target ship in order to hit it. That is, when you sweep your gaze to the left
beginning at the target ship (from own ship's taffrail), you see (1) the target ship's red flag, (2) the
target ship's blue flag, and (3) the point to aim at, or point "C," a considerable apparent distance
behind the target ship on the open water.
Acknowledgement: I used information from the "Hyper Physics Trajectories" web page by R. Nave, web site address
is http://hyperphysics.phy-astr.gsu.edu/hbase/traj.html although the formulae on that page are in the public
domain.
My email address is idouglas1304@yahoo.ca .
Program ARCSIN1.BAS. Microsoft GW-BASIC has no arc sine command or statement included, so I had to write one
of my own in order to do gun angle of elevation calculations. I obtained the algorithm for this from "Van
Nostrand's Scientific Encyclopedia," 6th edition, Considine et al editors, published 1983, page 1646a.
10 PRINT "Calculate ELANG = arc sin (RHEXPR)"
20 PRINT "ELANG means gun's angle of elevation; RHEXPR is input argument for the series approximation polynomial
used on right hand side of ELANG = arc sin (RHEXPR)"
30 PRINT:INPUT"Input argument for which arc sine is to be calculated ";RHEXPR
40 IF RHEXP>1 THEN PRINT "Input argument is out of range. ":GOTO 200
50 ELANG=RHEXPR+((RHEXPR^3)/6)+(.075*(RHEXPR^5))
60 ELANG=ELANG+(4.464285714285714D-02*(RHEXPR^7))
70 ELANG=ELANG+(3.038194444444444D-02*(RHEXPR^9))
80 ELANG=ELANG+(2.237215909090909D-02*(RHEXPR^11))
90 ELANG=ELANG+(1.735276442307692D-02*(RHEXPR^13))
100 ELANG=ELANG+(.01396484375#*(RHEXPR^15))
110 ELANG=ELANG+(1.155180089613971D-02*(RHEXPR^17))
120 ELANG=ELANG+(9.761609529194079D-03*(RHEXPR^19))
130 ELANG=ELANG+(8.390335809616816D-03*(RHEXPR^21))
140 ELANG=ELANG+(7.312525877035907D-03*(RHEXPR^23))
150 ELANG=ELANG+(6.447210311379962D-03*(RHEXPR^25))
160 ELANG=ELANG+(5.74003766956053D-03*(RHEXPR^27))
170 ELANG=ELANG+(.00515330968258#*(RHEXPR^29))
180 PRINT:PRINT "Arc sin of input argument (for this purpose, gun's angle of elevation):"
190 PRINT ELANG
200 PRINT:END
210 REM Did a loopback test in September 2009 and this program gave a reasonable
220 REM final answer given that there were only 7 decimal places in each answer
230 REM given by this program ARCSIN1.BAS.
240 REM I.e. did arcsin (0.523598775 radians) with this program;
250 REM got 0.5510696; did sin of this in GWBASIC using a direct statement;
260 REM got 0.5235989 radians via ARCSIN1.bas as the final answer.
Program TICTAC1.BAS
10 REM Begin flowchart box 1
20 CLEAR:CLS:PRINT "TICTAC program, that is, Time Interval"
30 PRINT "Compensating Tachymetry program, or"
40 PRINT "program to take data from two observations"
50 PRINT "of a moving target ship taken from"
60 PRINT "a moving ship and use"
70 PRINT "the data to determine the target ship's"
80 PRINT "starting and end X Y coordinates, azimuthal course"
90 PRINT "bearing, and speed. This program"
100 PRINT "is for IBM PC or similar computers"
110 PRINT "which do trigonometric calculations in radians."
120 PRINT "X y coordinates of own ship at start"
130 PRINT "are always assumed to be 50 nautical miles"
140 PRINT "by 50 nautical miles. Own ship's"
150 PRINT "true direction always assumed to be a"
160 PRINT "bearing of 0 degrees or 360 degrees."
170 PRINT "by Ian Douglas, year 2009."
180 PRINT "This program works for any azimuthal"
190 PRINT "bearing of target ship both at beginning and"
200 PRINT "end of time interval between the two observations."
210 REM End flowchart box 1
220 REM Begin flowchart box 2
230 PRINT "Enter the azimuthal bearing off own"
240 PRINT "ship of the target ship for first"
250 INPUT "observation (decimal degrees) :";ABT1:ABT1=ABT1/57.29577951#
260 PRINT "Enter the range to the target ship"
270 INPUT "for the first observation (feet): ";RT1
280 PRINT "Enter the azimuthal bearing off own"
290 PRINT "ship of the target ship for second"
300 INPUT "observation (decimal degrees) :";ABT2:ABT2=ABT2/57.29577951#
310 PRINT "Enter the range to the target ship"
320 INPUT "for the second observation (feet): ";RT2
330 PRINT "Enter the time elapsed in seconds between"
340 INPUT "The two observations ";TEBO
350 PRINT "Enter own ship's speed in"
360 INPUT "feet per second ";OSS
370 J1X=303805!:J1Y=303805!:REM Own
380 REM ship's starting coordinates in feet, i.e.
390 REM 50 nautical miles by 50 nautical miles
400 REM End flowchart box 2
410 REM Begin flowchart box 3
420 DOST=OSS*TEBO:REM DOST is distance
430 REM own ship travelled between the two observations
440 REM End flowchart box 3
450 REM Begin flowchart box 4
460 REM Calculate the coordinates of point J2, i.e.
470 REM end coordinates of own ship. J2's x
480 REM coordinate is simply J1's x coordinate
490 REM (303805 feet). J2's y coordinate is
500 REM J1's y coordinate (303805 feet) plus DOST.
510 J2X=J1X:J2Y=J1Y+DOST
520 REM End flowchart box 4
530 REM Begin flowchart box 5
540 REM Do a coordinated inverse (convert a vector and starting point
550 REM to an end point) to get the coordinates
560 REM of point K1, using point J1's coordinates
570 REM as the starting coordinates
580 REM Calculate end coordinates of vector
590 REM from angle, distance, and starting coordinates.
600 REM Distance is range to target
610 REM from first observation; azimuthal
620 REM bearing is that from first observation.
630 K1X=J1X+(RT1*SIN(ABT1)):K1Y=J1Y+(RT1*COS(ABT1))
640 REM End flowchart box 5
650 REM Begin flowchart box 6
660 REM Do another coordinate inverse (as before,
670 REM convert a vector and starting point to an end point)
680 REM to get the coordinates of point K2,
690 REM using the coordinates of point J2
700 REM as the starting coordinates and
710 REM observational data (range, bearing)
720 REM from point J2.
730 REM Calculate end coordinates of vector
740 REM from angle, distance, and starting coordinates.
750 REM Distance is range to target
760 REM from second observation; azimuthal
770 REM bearing is that from second observation.
780 K2X=J2X+(RT2*SIN(ABT2)):K2Y=J2Y+(RT2*COS(ABT2))
790 REM End flowchart box 6
800 REM Begin flowchart box 7
810 REM Do a coordinate inverse ("Convert 2
820 REM Cartesian points to vector") to get the
830 REM bearing and distance from K1 to K2.
840 DXAA=K2X-K1X:DYAA=K2Y-K1Y
850 DTAR=SQR((DXAA^2)+(DYAA^2)):REM Distance travelled by target
860 REM between the two observations.
870 IF DXAA=0 THEN GOTO 930
880 REM Above line prevents division by zero in next line
890 SLOPE=DYAA/DXAA:REM Slope
900 ANG1=ATN(SLOPE):REM Angle referred to standard position, calculated from slope
910 IF DXAA<0 THEN ANG1=ANG1+3.141592654#
920 IF DXAA>0 AND DYAA<0 THEN ANG1=ANG1+6.283185307#
930 IF DXAA=0 AND DYAA=0 THEN ANG1=0
940 IF DXAA=0 AND DYAA>0 THEN ANG1=1.570796327#
950 IF DXAA<0 AND DYAA=0 THEN ANG1=3.141592654#
960 IF DXAA=0 AND DYAA<0 THEN ANG1=4.712388981#
970 IF DXAA>0 AND DYAA=0 THEN ANG1=6.283185307#
980 AN1Z=6.283185308#-ANG1+1.570796327#
990 IF AN1Z>6.283185308# THEN AN1Z=AN1Z-6.283185308#
1000 REM Previous calculation gives azimuth from
1010 REM K1 to K2; 360 degrees minus angle plus 90 degrees
1020 REM End flowchart box 7
1030 REM Begin flowchart box 8
1040 REM Calculate target ship's speed
1050 TSS=DTAR/TEBO
1060 REM Recall that DTAR is distance travelled
1070 REM by target ship over TEBO,
1080 REM which is time elapsed between observations.
1090 REM End flowchart box 8
1100 REM Begin flowchart box 9
1110 PRINT "Target ship's speed (feet/sec): ";TSS
1120 PRINT "Target ship's azimuthal course"
1130 PRINT "bearing: ";AN1Z*57.29577951#;" decimal degrees of arc"
1140 PRINT "Target ship's starting coordinates (feet): ";K1X,K1Y
1150 PRINT "Target ship's end coordinates (feet): ";K2X,K2Y
1160 REM End flowchart box 9
1170 REM Begin flowchart box 10
1180 END
1190 REM End flowchart box 10
Program BATCON1.BAS
10 REM Begin flowchart box 1
20 CLEAR:CLS:PRINT "Battleship fire control program"
30 PRINT "for IBM PC or similar computers"
40 PRINT "which do trigonometric calculations in radians"
50 PRINT "by Ian Douglas, year 2007"
60 PRINT "This program works for any"
70 PRINT "azimuthal bearing of target ship"
80 PRINT "both at beginning and end of"
90 PRINT "shell's time of flight."
100 PRINT "Also works for any"
110 PRINT "bearing of point to aim at,"
120 PRINT "i.e. point 'C.'"
130 PRINT "X y coordinates of own ship at start"
140 PRINT "of problem are 50 nautical miles by 50 nautical miles."
150 PRINT "own ship's true direction always assumed to"
160 PRINT "be a bearing of 0 degrees or 360 degrees."
170 REM End flowchart box 1
180 REM Begin flowchart box 2
190 INPUT "Muzzle velocity in feet/sec ";MV
200 INPUT "Target ship's present range in feet ";TSR
210 INPUT "Target ship's true speed in feet/sec ";TSS
220 PRINT "Target ship's modified true plot direction of travel as an azimuthal bearing in"
230 INPUT "decimal degrees of arc "; TSD:TSD=TSD/57.29577951#
240 INPUT "Own ship's true speed in feet/sec";OSS
250 PRINT "Target ship's bearing off own ship"
260 INPUT "at time of gun's firing in decimal degrees of arc ";TSB:TSB=TSB/57.29577951#
270 G=32:REM acceleration due to gravity in feet per second per second
280 MAXRAN=((MV*MV)*SIN(1.570796327#))/G:REM Maximum theoretical range of gun
290 MAXRAN=MAXRAN*.85:REM Must use 85% (or thereabouts) of maximum theoretical range of gun
300 REM as maximum program-workable range of gun in order to prevent
310 REM possible looping, overflow, etc. errors associated with attempting
320 REM to fire gun at excessive ranges (albeit in some cases still under
330 REM maximum theoretical range of gun).
340 REM End flowchart box 2
350 REM Begin flowchart box 3
360 REM Calculate target ship's starting (i.e. initially sighted) coordinates in (feet, feet)
370 TSPX=TSR*SIN(TSB)+303805!
380 TSPY=TSR*COS(TSB)+303805!
390 REM Above addition of 303805 is number of feet in 50 nautical miles;
400 REM recall that own ship's position is (50 naut. miles, 50 naut. miles)
410 REM Calculate first approximation of gun's angle of elevation
420 THETA1=(TSR*G)/(MV*MV)
430 REM The following arc sine routine 'refines' THETA1 to an unnecessarily high
440 REM level of precision since all we need is an initial approximation
450 REM of the angle of the gun's elevation, and since we use a different
460 REM formula, later on, to calculate this angle; also the different formula
470 REM formula we use later on does not require calculation of y = arc sin (x).
480 THETA1=THETA1+((THETA1^3)/6)+(.075*(THETA1^5))
490 THETA1=THETA1+(4.464285714285714D-02*(THETA1^7))
500 THETA1=THETA1+(3.038194444444444D-02*(THETA1^9))
510 THETA1=THETA1+(2.237215909090909D-02*(THETA1^11))
520 THETA1=THETA1+(1.735276442307692D-02*(THETA1^13))
530 THETA1=THETA1+(.01396484375#*(THETA1^15))
540 THETA1=THETA1+(1.155180089613971D-02*(THETA1^17))
550 THETA1=THETA1+(9.761609529194079D-03*(THETA1^19))
560 THETA1=THETA1+(8.390335809616816D-03*(THETA1^21))
570 THETA1=THETA1+(7.312525877035907D-03*(THETA1^23))
580 THETA1=THETA1+(6.447210311379962D-03*(THETA1^25))
590 THETA1=THETA1+(5.74003766956053D-03*(THETA1^27))
600 THETA1=THETA1+(.00515330968258#*(THETA1^29))
610 THETA1=THETA1/2
620 REM Angle of elevation is theta1; used formula on p. 12 of notes.
630 PRINT "Target ship's starting (i.e. initially sighted) position in (feet, feet): "
640 PRINT TSPX,TSPY
650 PRINT "First approximation of gun's angle of elevation"
660 PRINT "(decimal degrees): ";THETA1*57.29577951#
670 REM Previous 3 lines are for debugging and error isolation only; delete or rem them out later
680 REM End flowchart box 3
690 REM Begin flowchart box 4
700 TOF1=2*((MV*SIN(THETA1))/G)
710 REM TOF1 is first approximation of shell's time of flight
720 PRINT "First approximation of shell's time of flight in seconds:"
730 PRINT TOF1
740 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
750 REM End flowchart box 4
760 REM Begin flowchart box 5
770 BEEP:REM Gets computer to beep to show that the program is looping back properly
780 REM on first (presumably wrong) guess of time of flight; need to use
790 REM BEEP since I cannot get any of my computers to 'listen' to the
800 REM scroll lock key on the keyboard(s). BEEP statement also delays operation
810 REM of program so you can see on screen whether proper reiterative looping
820 REM to get progressively better answers is taking place.
830 REM First step in this flowchart box is to get distance target ship
840 REM travelled during TOF (time of flight) calculated above.
850 DTST=TSS*TOF1:REM DTST is distance target ship travels during TOF1
860 REM (time of flight) calculated above
870 REM Now calculate how far own ship travelled during time of shell's flight
880 DOST=OSS*TOF1:REM DOST is how far own ship travelled during shell's
890 REM time of flight
900 REM Next step is to use section B of coordinate inverse program pp. 75-76
910 REM to calculate the end coordinates of the target ship at end of time
920 REM of flight, using azimuth of target ship, distance travelled, and
930 REM target's starting coordinates.
940 REM Input is target ship's initial sighted coordinates (TSPX, TSPY),
950 REM target ship's modified true plot direction of travel azimuthal bearing
960 REM TSD and distance target ship travels during shell's time of flight DTST.
970 REM Output is target ship's endpoint coordinates (TSEPX, TSEPY).
980 REM Calculate end coordinates
990 REM of vector, from angle, distance
1000 REM and starting coordinates
1010 TSEPX=TSPX+(DTST*SIN(TSD)):TSEPY=TSPY+(DTST*COS(TSD))
1020 REM End calculation of target ship's end coordinates
1030 REM Begin calculation of coordinates of point "C" (PCX, PCY)
1040 REM X coordinate of point "C" is simply x coordinate of target ship's end point
1050 REM at end of shell's time of flight, i.e. PCX=TSEPX
1060 PCX=TSEPX
1070 REM Y coordinate of point "C" is the y coordinate of target ship's endpoint at end of time of flight of shell TSEPY
1080 REM minus distance own ship travels during shell's time of flight
1090 PCY=TSEPY-DOST
1100 REM End calculation of coordinates of point "C"
1110 PRINT "Distance target ship travels during shell's time of flight (feet):"
1120 PRINT DTST
1130 PRINT "Distance own ship travels during shell's time of flight (feet):"
1140 PRINT DOST
1150 PRINT "Target ship's end coordinates (feet, feet):"
1160 PRINT TSEPX,TSEPY
1170 PRINT "Coordinates of point 'C' (feet, feet):"
1180 PRINT PCX,PCY
1190 REM Previous 8 lines are for debugging and error isolation only; delete or rem them out later
1200 REM End flowchart box 5
1210 REM Begin flowchart box 6
1220 REM Calculate straight line distance between own ship at time of gun's firing
1230 REM and point "C," DOSPC
1240 REM Recall that own ship's coordinates at time of gun's firing are
1250 REM (303805 feet, 303805 feet) (i.e. 50 nautical miles by 50 nautical miles)
1260 DOSPC=SQR(((PCX-303805!)^2)+((PCY-303805!)^2))
1270 IF DOSPC>MAXRAN THEN PRINT "Target ship is out of range. ":GOTO 1760
1280 PRINT "Distance between own ship and point 'C' when gun is fired (feet):"
1290 PRINT DOSPC
1300 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
1310 REM End flowchart box 6
1320 REM Begin flowchart box 7
1330 THETA2=(DOSPC*G)/(MV*MV)
1340 TOF2=2*(MV*SIN(THETA2))/G
1350 PRINT "New elevation and time of flight (decimal degrees of arc and seconds of time):"
1360 PRINT THETA2*57.29577951#,TOF2
1370 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
1380 REM End flowchart box 7
1390 REM Begin flowchart box 8
1400 IF (ABS(TOF2-TOF1))<.001 THEN GOTO 1490
1410 REM Note that the above line of code gives 1/1000 of a second accuracy!! I have found through
1420 REM experimentation that seeking higher accuracy will in some cases cause
1430 REM infinite looping of the program in a search for unnecessarily high accuracy.
1440 REM It is debatable whether in this program you need even 1/100 of a second accuracy.
1450 REM End flowchart box 8
1460 REM Begin flowchart box 9
1470 TOF1=TOF2:GOTO 760
1480 REM End flowchart box 9
1490 REM Begin flowchart box 10
1500 REM End flowchart box 10
1510 REM Begin flowchart box 11
1520 REM Calculate gun's bearing using coordinates of point "C"
1530 REM From flowchart box 5, and using coordinates of own ship at time of gun's firing;
1540 REM Convert points to vector
1550 REM 1st point is coordinates of own ship at gun's time of firing (303805 feet, 303805 feet)
1560 REM 2nd point is coordinates of point "C," (PCX, PCY)
1570 DXAA=PCX-303805!:DYAA=PCY-303805!
1580 IF DXAA=0 THEN GOTO 1640
1590 REM Above line prevents division by zero in next line
1600 SLOPE=DYAA/DXAA:REM Slope
1610 ANGLE#=ATN(SLOPE):REM Angle referred to standard position, calculated from slope
1620 IF DXAA<0 THEN ANGLE#=ANGLE#+3.141592654#
1630 IF DXAA>0 AND DYAA<0 THEN ANGLE#=ANGLE#+6.283185307#
1640 IF DXAA=0 AND DYAA=0 THEN ANGLE#=0
1650 IF DXAA=0 AND DYAA>0 THEN ANGLE#=1.570796327#
1660 IF DXAA<0 AND DYAA=0 THEN ANGLE#=3.141592654#
1670 IF DXAA=0 AND DYAA<0 THEN ANGLE#=4.712388981#
1680 IF DXAA>0 AND DYAA=0 THEN ANGLE#=6.283185307#
1690 REM End flowchart box 11
1700 REM Begin flowchart box 12
1710 REM Gun's angle of elevation is previously calculated THETA2
1720 PRINT "Gun's azimuth to shoot along: ";GUNAZ#*57.29577951#;" decimal degrees of arc"
1730 PRINT "Gun's angle of elevation: ";THETA2*57.29577951#;" decimal degrees of arc"
1740 REM If latest value for time of flight is to be printed here, that value is TOF2.
1750 REM End flowchart box 12
1760 REM Begin flowchart box 13
1770 END
1780 REM End flowchart box 13
Program UPART.BAS: this program uses a trial and error method to calculate the angle of elevation of a gun needed
to hit a target of higher elevation than the gun. The formula to do this is a transcendental function (similar
to Kepler's Equation) and cannot be solved directly. This program assumes that the target is within range and will
simply crash or go into an infinite loop if the target whose coordinates are entered into the program is out of
range. See the Appendix for a discussion of how to determine whether a target of higher elevation than the gun
is within the gun's range, and for a more advanced version of this program which traps such out-of-range-of-the-gun
situations. Also note that this program goes into an infinite loop if the fire and hit solution involves the gun
having an angle of elevation of greater than about 78 degrees (the true value is higher but I have not determined
this exactly so to be safe I have given a lower and rounded off figure). For example, if your gun has a muzzle velocity
of 611.8823 feet per second and the target is separated from the gun vertically and horizontally by 4900 feet and
1460 feet, then the program will go into an infinite loop; but if the horizontal separation increases
to 1470 feet then the program will tell you to elevate the gun to 78.4921 decimal degrees above horizontal. In other
words, the program will not work for targets directly overhead and above the gun or nearly so. I hope at some time
to address this great weakness but for the time being the program's end user must beware of this flaw. On the other
hand, note that this weakness creeps in when the gun is firing only 12 degrees of arc--or less--from vertical.
I am still debating how to handle the last mentioned types of situation; such instances mean that the program
has a blind spot where a fire and hit solution is physically possible but not practically computable. There are
two ways to handle such instances. One is to get the program to end if it tries too many times to find an answer
to the question of what the gun's angle of elevation should be to hit the target. The other method is to get the
program to end if it starts to use, as trial angles of elevation, angles somewhere between 75 and 78 degrees. The
first method is virtually guaranteed to work but might take some time to be processed by a program; the second method
is potentially faster but could be dangerous, as the program might accidentally cross the threshold of "too high" an
angle of elevation by too great an increment. This having occurred, the program would run unchecked into the domain
where it cannot use practical means to find an answer--plunging itself into an infinite loop. At some point in the
future I hope to make a decision and incorporate the decision into a revised program UPART.BAS and other programs on
this web site related to it, e.g. AIRBAT.BAS.
10 REM Begin flowchart box E1 (begin/end marker)
20 REM End flowchart box E1; the flowchart numbers in this program correspond to flowchart on p. 251 of notes.
30 REM Begin flowchart box E2
40 PRINT "Program to use a variant of successive approximation to get angle of elevation"
50 PRINT "of a gun firing at a target separated from the gun significantly"
60 PRINT "separated from it both vertically and horizontally (e.g. an aircraft) using equation"
70 PRINT "211a on page 211 of notes; by Ian Douglas, year 2010. This program also works"
80 PRINT "for targets below gun. This program is for upward-only trajectories if target is above gun"
90 PRINT "and not for upward-to-peak-to-downward trajectories."
100 REM Program is intended to make the problem manageable by taking the maximum possible shell height the gun
110 REM can attain and obtaining a constant to multiply 90 degrees of arc (1.570796327 radians) (gun's maximum angle of
120 REM elevation) by to get any given possible height the shell can attain; this enables a resultant vertical
130 REM separation of the target from gun to be compared to an actual vertical separation attained by a trial angle
140 REM of elevation of gun, for comparison purposes; also this enables the successive approximation process to occur.
150 REM For this reason, I do not recommend using this program to compute data concerning guns with greater
160 REM than 10,000 feet per second muzzle velocities; but even the most powerful guns have about 1/4 that.
170 REM Also note that the program will go into an infinite loop if the target is directly above gun, even though
180 REM a fire and hit solution is possible in such a case.
190 PRINT "Vertical separation (Y) of target from gun (feet) "
200 PRINT "(negative for a target below gun)"
210 INPUT VESEP
220 INPUT "Horizontal separation (in XY plane) of target from gun (feet) ";HOSEP
230 INPUT "Gun's muzzle velocity (feet per second) ";MV
240 G=32:G1=G/2:REM Acceleration due to gravity in feet per second per second; also half this for use in
250 REM instances of equation 211a.
260 REM End flowchart box E2
270 REM Begin flowchart box E3
280 REM Calculate the shell's maximum attainable height.
290 MAXHT=(MV*MV)/(2*G)
300 REM Calculate the constant mentioned in remarks for flowchart box E2.
310 CONST=MAXHT/1.570796327#
320 REM Calculate the minimum angle of elevation of gun needed to reach HOSEP.
330 MIELAN=(HOSEP*G)/(MV*MV)
340 REM The following is an arc sine routine with the final angle divided by 2 for the purpose of the
350 REM above formula to get the minimum necessary angle of elevation given HOSEP, muzzle velocity, and G.
360 MIELAN=MIELAN+((MIELAN^3)/6)+(.075*(MIELAN^5))
370 MIELAN=MIELAN+(4.464285714285714D-02*(MIELAN^7))
380 MIELAN=MIELAN+(3.038194444444444D-02*(MIELAN^9))
390 MIELAN=MIELAN+(2.237215909090909D-02*(MIELAN^11))
400 MIELAN=MIELAN+(1.735276442307692D-02*(MIELAN^13))
410 MIELAN=MIELAN+(.01396484375#*(MIELAN^15))
420 MIELAN=MIELAN+(1.155180089613971D-02*(MIELAN^17))
430 MIELAN=MIELAN+(9.761609529194079D-03*(MIELAN^19))
440 MIELAN=MIELAN+(8.390335809616816D-03*(MIELAN^21))
450 MIELAN=MIELAN+(7.312525877035907D-03*(MIELAN^23))
460 MIELAN=MIELAN+(6.447210311379962D-03*(MIELAN^25))
470 MIELAN=MIELAN+(5.74003766956053D-03*(MIELAN^27))
480 MIELAN=MIELAN+(.00515330968258#*(MIELAN^29))
490 MIELAN=MIELAN/2
500 REM Angle of elevation is MIELAN; used formula on p. 12 of notes.
510 REM Set shell's minimum possible height as zero feet; set gun's maximum angle of elevation as 90 degrees of arc.
520 REM End flowchart box E3
530 REM Begin flowchart box E4
540 REM Set first trial angle of elevation of gun as (45 degrees of arc minus minimum angle of elevation of gun
550 REM to reach HOSEP)/2; the reason I use 45 and not 90 degrees here is that as soon as you use an angle greater
560 REM than 45 degrees you meet up with the possibility of the shell falling short of HOSEP.
570 ELANG=(.785398163#-MIELAN)/2
580 REM End flowchart box E4
590 REM Begin flowchart box E5
600 REM Calculate trial vertical (Y) separation resulting from ELANG and inputted values from flowchart box E2
610 REM (except VESEP) using equation 211a; calculate difference between trial vertical separation and VESEP.
620 A1=MV*SIN(ELANG):A2=MV*COS(ELANG)
630 TRIZ=(HOSEP*A1)/A2
640 TRIZ=TRIZ-G1*((HOSEP*HOSEP)/(A2*A2))
650 REM PRINT "Trial Z and trial gun's angle of elevation are ";TRIZ,ELANG*57.29577951#:REM For debugging only
660 REM End flowchart box E5
670 REM Begin flowchart box E6
680 IF ABS(TRIZ-VESEP)<.01 THEN GOTO 770:REM That is, go to flowchart box E8. Note that this represents a 1/100
690 REM of a foot accuracy of the shell's striking point in relation to the program's perception of the target,
700 REM where the target is an aircraft! Compare this with the size of any conventional aircraft.
710 REM I have tried to get the program to achieve a higher level of precision, but this only caused
720 REM endless looping due to the problem of attainable accuracy of TRIZ and ELANG, loss of digits, etc.
730 REM End flowchart box E6
740 REM Begin flowchart box E7
750 ELANG=ELANG-(((TRIZ-VESEP)/2)/CONST):GOTO 590:REM That is, go to flowchart box E5
760 REM End flowchart box E7
770 REM Begin flowchart box E8
780 PRINT:PRINT "Gun's angle of elevation from horizontal (depression below horizontal negative) = ";
790 PRINT ELANG*57.29577951#;" decimal degrees of arc":PRINT
800 REM End flowchart box E8
810 REM Begin flowchart box E9 (begin/end marker)
820 END
830 REM End flowchart box E9
Program ATICTAC.BAS
10 REM Begin flowchart box 1
20 CLEAR:CLS:PRINT "ATICTAC program, that is, Aircraft Time Interval"
30 PRINT "Compensating Tachymetry program, or"
40 PRINT "program to take data from two observations"
50 PRINT "of a moving target aircraft taken from"
60 PRINT "a moving ship and use"
70 PRINT "the data to determine the target aircraft's"
80 PRINT "starting and end XYZ coordinates, azimuthal course"
90 PRINT "bearing, speed, and rate of climb (descent negative). This program"
100 PRINT "is for IBM PC or similar computers"
110 PRINT "which do trigonometric calculations in radians."
120 REM Note that X and Y coordinates are
130 REM ground coordinates and Z
140 REM is altitiude above own ship.
150 PRINT:PRINT "X Y coordinates of own ship at start"
160 PRINT "are always assumed to be 50 nautical miles"
170 PRINT "by 50 nautical miles. Own ship's"
180 PRINT "true direction always assumed to be a"
190 PRINT "bearing of 0 degrees or 360 degrees."
200 PRINT:PRINT "by Ian Douglas, year 2010."
210 PRINT:PRINT "This program works for any azimuthal"
220 PRINT "bearing of target aircraft both at beginning and"
230 PRINT "end of time interval between the two observations."
240 REM This program assumes altitude of own ship always to be 0 feet and works only for positive altitudes of aircraft.
250 REM End flowchart box 1
260 REM Begin flowchart box 2
270 PRINT:PRINT "Enter the azimuthal bearing off own"
280 PRINT "ship of the target aircraft for first"
290 INPUT "observation (decimal degrees) :";ABT1:ABT1=ABT1/57.29577951#
300 PRINT "Enter the straight line of sight range to the target aircraft"
310 INPUT "for the first observation (feet): ";RT1
320 PRINT "Enter the observed angle of elevation of the target aircraft in decimal degrees"
330 INPUT "above horizontal for the first observation";ELANG1:ELANG1=ELANG1/57.29577951#
340 PRINT:PRINT "Enter the azimuthal bearing off own"
350 PRINT "ship of the target aircraft for second"
360 INPUT "observation (decimal degrees) :";ABT2:ABT2=ABT2/57.29577951#
370 PRINT "Enter the straight line of sight range to the target aircraft"
380 INPUT "for the second observation (feet): ";RT2
390 PRINT "Enter the observed angle of elevation of the target aircraft in decimal degrees"
400 INPUT "above horizontal for the second observation";ELANG2:ELANG2=ELANG2/57.29577951#
410 PRINT:PRINT "Enter the time elapsed in seconds between"
420 INPUT "The two observations ";TEBO
430 PRINT:PRINT "Enter own ship's speed in"
440 INPUT "feet per second ";OSS
450 J1X=303805!:J1Y=303805!:REM Own
460 REM Own ship's starting coordinates in feet, i.e.
470 REM 50 nautical miles by 50 nautical miles
480 REM End flowchart box 2
490 REM Begin flowchart box 3
500 ALT1=RT1*(SIN(ELANG1)):ALT2=RT2*(SIN(ELANG2)):REM Computes altitudes of aircraft for both observations
510 RCLMB=(ALT2-ALT1)/TEBO:REM rate of climb of aircraft (descent negative)
520 HRAN1=RT1*(COS(ELANG1)):HRAN2=RT2*(COS(ELANG2)):REM Computes horizontal ground track ranges to target aircraft for both observations
530 DOST=OSS*TEBO:REM DOST is distance
540 REM own ship travelled between the two observations
550 REM End flowchart box 3
560 REM Begin flowchart box 4
570 REM Calculate the coordinates of point J2, i.e.
580 REM end coordinates of own ship. J2's x
590 REM coordinate is simply J1's x coordinate
600 REM (303805 feet). J2's y coordinate is
610 REM J1's y coordinate (303805 feet) plus DOST.
620 J2X=J1X:J2Y=J1Y+DOST
630 REM End flowchart box 4
640 REM Begin flowchart box 5
650 REM Do a coordinate inverse (convert a vector and starting point
660 REM to an end point) to get the coordinates
670 REM of point K1, using point J1's coordinates
680 REM as the starting coordinates
690 REM Calculate end coordinates of vector
700 REM from angle, distance, and starting coordinates.
710 REM Distance is horizontal ground track range to target aircraft
720 REM from first observation; azimuthal
730 REM bearing is that from first observation.
740 K1X=J1X+(HRAN1*SIN(ABT1)):K1Y=J1Y+(HRAN1*COS(ABT1)):K1Z=ALT1
750 REM End flowchart box 5
760 REM Begin flowchart box 6
770 REM Do another coordinate inverse (as before,
780 REM convert a vector and a starting point to an end point)
790 REM to get the coordinates of point K2,
800 REM using the coordinates of point J2
810 REM as the starting coordinates and
820 REM observational data (range, bearing)
830 REM from point J2.
840 REM Calculate end coordinates of vector
850 REM from angle, distance, and starting coordinates.
860 REM Distance is horizontal ground track range to target
870 REM from second observation; azimuthal
880 REM bearing is that from second observation.
890 K2X=J2X+(HRAN2*SIN(ABT2)):K2Y=J2Y+(HRAN2*COS(ABT2)):K2Z=ALT2
900 REM End flowchart box 6
910 REM Begin flowchart box 7
920 REM Do a coordinate inverse ("Convert 2
930 REM Cartesian points to vector") to get the
940 REM bearing and distance from K1 to K2.
950 DXAA=K2X-K1X:DYAA=K2Y-K1Y:DZAA=K2Z-K1Z
960 DTAR=SQR((DXAA^2)+(DYAA^2)):REM Ground track distance travelled by target aircraft
970 REM between the two observations.
980 IF DXAA=0 THEN GOTO 1040
990 REM Above line prevents division by zero in next line
1000 SLOPE=DYAA/DXAA:REM Slope
1010 ANG1=ATN(SLOPE):REM Angle referred to standard position, calculated from slope
1020 IF DXAA<0 THEN ANG1=ANG1+3.141592654#
1030 IF DXAA>0 AND DYAA<0 THEN ANG1=ANG1+6.283185307#
1040 IF DXAA=0 AND DYAA=0 THEN ANG1=0
1050 IF DXAA=0 AND DYAA>0 THEN ANG1=1.570796327#
1060 IF DXAA<0 AND DYAA=0 THEN ANG1=3.141592654#
1070 IF DXAA=0 AND DYAA<0 THEN ANG1=4.712388981#
1080 IF DXAA>0 AND DYAA=0 THEN ANG1=6.283185307#
1090 AN1Z=6.283185308#-ANG1+1.570796327#
1100 IF AN1Z>6.283185308# THEN AN1Z=AN1Z-6.283185308#
1110 REM Previous calculation gives azimuth from
1120 REM K1 to K2; 360 degrees minus angle plus 90 degrees
1130 REM End flowchart box 7
1140 REM Begin flowchart box 8
1150 REM Calculate target aircraft's ground track speed
1160 TAS=DTAR/TEBO
1170 REM Recall that DTAR is ground track distance travelled
1180 REM by target aircraft over TEBO,
1190 REM which is time elapsed between observations.
1200 REM End flowchart box 8
1210 REM Begin flowchart box 9
1220 PRINT:PRINT "Target aircraft's ground track speed (feet/sec): ";TAS
1230 PRINT "Target aircraft's azimuthal course"
1240 PRINT "bearing: ";AN1Z*57.29577951#;" decimal degrees of arc"
1250 PRINT "Target aircraft's rate of"
1260 PRINT "climb (descent negative) (feet per second) ";RCLMB
1270 PRINT "Target aircraft's starting XYZ coordinates (feet): ";K1X,K1Y,K1Z
1280 PRINT "Target aircraft's end XYZ coordinates (feet): ";K2X,K2Y,K2Z
1290 REM End flowchart box 9
1300 REM Begin flowchart box 10
1310 END
1320 REM End flowchart box 10
Program AIRBAT.BAS : This program assumes that the target is within range both at the outset of the problem
(initial acqisition of the target) and at the end of the problem (hitting the target). The program will
simply crash or go into an infinite loop if the target whose coordinates are entered into the program is out of
range. See the Appendix for a discussion of how to determine whether a target of higher elevation than the gun
is within the gun's range. Also note that this program has a weakness also present in UPART.BAS, which is that it
goes into an infinite loop if the fire and hit solution involves the gun having an angle of elevation of greater than
about 78 degrees: see the introductory notes for the program UPART.BAS above for a fuller explanation.
10 REM Begin flowchart box F1
20 CLEAR:CLS:PRINT "Antiaircraft gun fire control program for firing at a moving aircraft from a moving ship."
30 PRINT "This program is for IBM PC or similar computers"
40 PRINT "which do trigonometric calculations in radians"
50 PRINT "by Ian Douglas, year 2010"
60 PRINT "This program works for any azimuthal bearing of target aircraft"
70 PRINT "both at beginning and end of shell's time of flight."
80 PRINT "Also works for any bearing of point to aim at, i.e. point 'C.'"
90 PRINT "Note that this program will not point the gun directly at point 'C' but to a point directly above it,"
100 PRINT "not named or calculated in the program, and let the shell fall to point 'C' by gravity."
110 PRINT "X Y Z coordinates of own ship at start"
120 PRINT "of problem are 50 nautical miles by 50 nautical miles by 0."
130 PRINT "Own ship's true direction always assumed to"
140 PRINT "be a bearing of 0 degrees or 360 degrees."
150 REM This program assumes that it is working at or near the earth's surface, where the shell flies at the target
160 REM under gravitational influences.
170 PRINT "This program works for targets whose starting and/or ending at elevations are above or below own ship;"
180 PRINT "likewise for the shell's point of impact. If gun must be depressed to hit target, then gun's"
190 PRINT "angle of elevation is negative."
200 REM End flowchart box F1
210 REM Begin flowchart box F2
220 INPUT "Muzzle velocity in feet/sec ";MV
230 INPUT "Target aircraft's present straight line-of-sight range in feet ";TAR
240 INPUT "Target aircraft's true 'ground track' speed in feet/sec ";TAS
250 PRINT "Target aircraft's modified true plot 'ground track' direction of travel as an"
260 INPUT "azimuthal bearing in decimal degrees of arc "; TAD:TAD=TAD/57.29577951#
270 PRINT "Target aircraft's first observed angle of elevation above horizontal"
280 INPUT "in decimal degrees of arc ";TAAE:TAAE=TAAE/57.29577951#
290 PRINT "Target aircraft's rate of climb in feet/sec"
300 INPUT "(descent negative) ";TARC
310 INPUT "Own ship's true speed in feet/sec";OSS
320 PRINT "Target aircraft's bearing off own ship"
330 INPUT "at time of gun's firing in decimal degrees of arc ";TAB:TAB=TAB/57.29577951#
340 G=32:G1=G/2:REM Acceleration due to gravity in feet per second per second; also half this for use in
350 REM instances of equation 211a.
360 REM End flowchart box F2
370 REM Begin flowchart box F3
380 REM Calculate the shell's maximum attainable height for the purposes of calculating the value of the
390 REM constant mentioned in remarks for flowchart box E2 for subroutine 'A'; also do that second calculation.
400 REM I have moved the function of these calculations to here since subroutine 'A'
410 REM works inside a larger program, and since these calculations need only be done once.
420 MAXHT=(MV*MV)/(2*G):CONST=MAXHT/1.570796327#
430 REM Calculate target aircraft's flat range (straight line of sight range reduced to XY plane) and then
440 REM calculate beginning (i.e. initially sighted) coordinates in (feet, feet)
450 FLRANG=TAR*COS(TAAE)
460 TABX=(FLRANG*SIN(TAB))+303805!
470 TABY=(FLRANG*COS(TAB))+303805!
480 TABZ=TAR*SIN(TAAE)
490 REM Above addition of 303805 is number of feet in 50 nautical mile;
500 REM recall that own ship's position is (50 naut. mile, 50 naut. mile, 0)
510 PRINT "Target aircraft's starting (i.e. initially sighted) XYZ position in (feet, feet, feet): "
520 PRINT TABX,TABY,TABZ
530 REM Previous 3 lines are for debugging and error isolation only; delete or rem them out later
540 REM Calculate initial horizontal and vertical separation between own ship at time of gun's firing
550 REM and target aircraft (i.e. initial values of HOSEP and VESEP)
560 REM Recall that own ship's coordinates at time of gun's firing are
570 REM (303805 feet, 303805 feet, 0 feet) (i.e. 50 nautical miles by 50 nautical miles by 0)
580 HOSEP=SQR(((TABX-303805!)^2)+(TABY-303805!)^2)
590 VESEP=TABZ
600 PRINT "Horizontal and vertical distances between own ship and target aircraft:"
610 PRINT HOSEP,VESEP
620 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
630 GOSUB 1770:TOF1=TOF2:REM That is, go to subroutine A to get latest approximation of shell's time of flight TOF1, and
640 REM first approximation of gun's angle of elevation (this last not needed in first iterations of program,
650 REM but I wanted to have a single 'package' subroutine to do the necessary multiple tasks in later
660 REM iterations of the problem). Note that since the subroutine calls the latest time of flight TOF2, and
670 REM since we are in a position in the program where we need only the first approximation of time of flight
680 REM (and use this in the initial calculations immediately following this line of code) we name this result TOF1.
690 REM Later we will stick with the name the subroutine calls this time of flight, i.e. TOF2, unchanged.
700 PRINT "First approximation of shell's time of flight in seconds:"
710 PRINT TOF1
720 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
730 REM End flowchart box F3
740 REM Begin flowchart box F4 (off page connector)
750 REM End flowchart box F4 (off page connector)
760 REM Begin flowchart box F5
770 BEEP:REM Gets computer to beep to show that the program is looping back properly
780 REM on first (presumably wrong) guess of time of flight; need to use
790 REM BEEP since I cannot get any of my computers to 'listen' to the
800 REM scroll lock key on the keyboard(s). BEEP statement also delays operation
810 REM of program so you can see on screen whether proper reiterative looping
820 REM to get progressively better answers is taking place.
830 REM First step in this flowchart box is to get distance target aircraft
840 REM travelled during TOF (time of flight) calculated above.
850 DTAT=TAS*TOF1:REM DTAT is 'ground track' distance target aircraft travels during TOF1
860 REM (time of flight) calculated above
870 REM Now calculate how far own ship travelled during time of shell's flight
880 DOST=OSS*TOF1:REM DOST is 'ground track'distance own ship travelled during shell's
890 REM time of flight
900 REM Next step is to use section B of coordinate inverse program pp. 75-76
910 REM to calculate the end coordinates of the target aircraft at end of time
920 REM of flight, using azimuth of target aircraft, distance travelled, and
930 REM target aircraft's starting coordinates.
940 REM Input is target aircraft's initial sighted coordinates (TABX, TABY, TABZ),
950 REM target aircraft's modified true plot direction of travel azimuthal bearing
960 REM TAD and distance target aircraft travels during shell's time of flight DTAT.
970 REM Output is target aircraft's endpoint coordinates (TAEPX, TAEPY, TAEPZ).
980 REM Calculate end coordinates
990 REM of vector, from angle, distance
1000 REM and starting coordinates
1010 TAEPX=TABX+(DTAT*SIN(TAD)):TAEPY=TABY+(DTAT*COS(TAD)):TAEPZ=TABZ+(TARC*TOF1)
1020 REM End calculation of target aircraft's end coordinates
1030 REM Begin calculation of coordinates of point "C" (PCX, PCY, PCZ)
1040 REM X coordinate of point "C" is simply x coordinate of target aircraft's end point
1050 REM at end of shell's time of flight, i.e. PCX=TAEPX
1060 PCX=TAEPX
1070 REM Y coordinate of point "C" is the y coordinate of target aircraft's endpoint at end of time of flight of shell TAEPY
1080 REM minus distance own ship travels during shell's time of flight
1090 PCY=TAEPY-DOST
1100 REM Z coordinate of point "C" is simply target aircraft's altitude at end of shell's time of flight, i.e. TAEPZ
1110 PCZ=TAEPZ
1120 REM End calculation of coordinates of point "C"
1130 PRINT "'Ground track' distance target aircraft travels during shell's time of flight (feet):"
1140 PRINT DTAT
1150 PRINT "'Ground track' distance own ship travels during shell's time of flight (feet):"
1160 PRINT DOST
1170 PRINT "Target aircraft's XYZ end coordinates (feet, feet, feet):"
1180 PRINT TAEPX,TAEPY,TAEPZ
1190 PRINT "XYZ coordinates of point 'C' (feet, feet, feet):"
1200 PRINT PCX,PCY,PCZ
1210 REM Previous 8 lines are for debugging and error isolation only; delete or rem them out later
1220 REM End flowchart box F5
1230 REM Begin flowchart box F6
1240 REM Calculate horizontal and vertical separation between own ship at time of gun's firing
1250 REM and point "C," (i.e. HOSEP and VESEP)
1260 REM Recall that own ship's coordinates at time of gun's firing are
1270 REM (303805 feet, 303805 feet, 0 feet) (i.e. 50 nautical miles by 50 nautical miles by 0)
1280 HOSEP=SQR(((PCX-303805!)^2)+(PCY-303805!)^2)
1290 VESEP=PCZ
1300 PRINT "Horizontal and vertical distances between own ship and point 'C' when gun is fired (feet):"
1310 PRINT HOSEP,VESEP
1320 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
1330 GOSUB 1770:REM That is, go to subroutine A to get latest approximation of shell's time of flight TOF2, and
1340 REM latest approximation of gun's angle of elevation.
1350 REM End flowchart box F6
1360 REM Begin flowchart box F7
1370 IF (ABS(TOF2-TOF1))<.001 THEN GOTO 1470:REM That is, go to flowchart box F9
1380 REM Note that the above line of code gives 1/1000 of a second accuracy!! I have found through
1390 REM experimentation that seeking higher accuracy will in some cases cause
1400 REM infinite looping of the program in a search for unnecessarily high accuracy.
1410 REM It is debatable whether in this program you need even 1/100 of a second accuracy. To make the program get an
1420 REM answer in less time, reduce accuracy by ten or a hundred fold in the line of code 5 lines before this one.
1430 REM End flowchart box F7
1440 REM Begin flowchart box F8
1450 TOF1=TOF2:GOTO 760:REM That is, go to flowchart box F5
1460 REM End flowchart box F8
1470 REM Begin flowchart box F9
1480 REM Calculate gun's bearing (not angle of elevation yet) using only X Y coordinates of point "C"
1490 REM From flowchart box 4, and using coordinates of own ship at time of gun's firing;
1500 REM Convert points to vector
1510 REM 1st point is coordinates of own ship at gun's time of firing (303805 feet, 303805 feet)
1520 REM 2nd point is coordinates of point "C," (PCX, PCY)
1530 DXAA=PCX-303805!:DYAA=PCY-303805!
1540 IF DXAA=0 THEN GOTO 1600
1550 REM Above line prevents division by zero in next line
1560 SLOPE=DYAA/DXAA:REM Slope
1570 ANGLE#=ATN(SLOPE):REM Angle referred to standard position, calculated from slope
1580 IF DXAA<0 THEN ANGLE#=ANGLE#+3.141592654#
1590 IF DXAA>0 AND DYAA<0 THEN ANGLE#=ANGLE#+6.283185307#
1600 IF DXAA=0 AND DYAA=0 THEN ANGLE#=0
1610 IF DXAA=0 AND DYAA>0 THEN ANGLE#=1.570796327#
1620 IF DXAA<0 AND DYAA=0 THEN ANGLE#=3.141592654#
1630 IF DXAA=0 AND DYAA<0 THEN ANGLE#=4.712388981#
1640 IF DXAA>0 AND DYAA=0 THEN ANGLE#=6.283185307#
1650 GUNAZ#=6.283185307#-ANGLE#+1.570796327#:REM 360 degrees minus angle plus 90 degrees
1660 IF GUNAZ#>6.283185307# THEN GUNAZ#=GUNAZ#-6.283185307#
1670 REM End flowchart box F9
1680 REM Begin flowchart box F10
1690 REM Gun's angle of elevation is previously calculated ELANG
1700 PRINT "Gun's azimuth to shoot along: ";GUNAZ#*57.29577951#;" decimal degrees of arc"
1710 PRINT "Gun's angle of elevation: "; ELANG*57.29577951#;" decimal degrees of arc"
1720 REM If latest value for time of flight is to be printed here, that value is TOF2.
1730 REM End flowchart box F10
1740 REM Begin flowchart box F11
1750 END
1760 REM End flowchart box F11
1770 REM Begin flowchart box E1 (begin/end marker)
1780 REM End flowchart box E1; the flowchart numbers in this program correspond to flowchart on p. 251 of notes
1790 REM since this subroutine is simply (for the most part, except for the time of flight calculations
1800 REM at the end) a repeat of the program UPART.BAS.
1810 REM Begin flowchart box E2
1820 REM *** Subroutine 'A' ***, to use a variant of successive approximation to get angle of elevation
1830 REM of a gun firing at a target separated from the gun significantly
1840 REM separated from it both vertically and horizontally (e.g. an aircraft) using equation
1850 REM 211a on page 211 of notes; by Ian Douglas, year 2010. This program also works
1860 REM for targets below gun. This program is for upward-only trajectories if target is above gun
1870 REM and not for upward-to-peak-to-downward trajectories. This subroutine also calculates the
1880 REM shell's time of flight.
1890 REM Program is intended to make the problem manageable by taking the maximum possible shell height the gun
1900 REM can attain and obtaining a constant to multiply 90 degrees of arc (1.570796327 radians) (gun's maximum angle of
1910 REM elevation) by to get any given possible height the shell can attain; this enables a resultant vertical
1920 REM separation of the target from gun to be compared to an actual vertical separation attained by a trial angle
1930 REM of elevation of gun, for comparison purposes; also this enables the successive approximation process to occur.
1940 REM For this reason, I do not recommend using this program to compute data concerning guns with greater
1950 REM than 10,000 feet per second muzzle velocities; but even the most powerful guns have about 1/4 that.
1960 REM End flowchart box E2
1970 REM Begin flowchart box E3
1980 REM If this subroutine were working on its own, I would have put the calculations for the shell's maximum
1990 REM attainable height and for the constant needed in this subroutine here; but since this subroutine
2000 REM functions inside a larger program, and since these calculations need only be done once, I have
2010 REM left them to flowchart box F3.
2020 REM Calculate the minimum angle of elevation of gun needed to reach HOSEP.
2030 MIELAN=(HOSEP*G)/(MV*MV)
2040 REM The following is an arc sine routine with the final angle divided by 2 for the purpose of the
2050 REM above formula to get the minimum necessary angle of elevation given HOSEP, muzzle velocity, and G.
2060 MIELAN=MIELAN+((MIELAN^3)/6)+(.075*(MIELAN^5))
2070 MIELAN=MIELAN+(4.464285714285714D-02*(MIELAN^7))
2080 MIELAN=MIELAN+(3.038194444444444D-02*(MIELAN^9))
2090 MIELAN=MIELAN+(2.237215909090909D-02*(MIELAN^11))
2100 MIELAN=MIELAN+(1.735276442307692D-02*(MIELAN^13))
2110 MIELAN=MIELAN+(.01396484375#*(MIELAN^15))
2120 MIELAN=MIELAN+(1.155180089613971D-02*(MIELAN^17))
2130 MIELAN=MIELAN+(9.761609529194079D-03*(MIELAN^19))
2140 MIELAN=MIELAN+(8.390335809616816D-03*(MIELAN^21))
2150 MIELAN=MIELAN+(7.312525877035907D-03*(MIELAN^23))
2160 MIELAN=MIELAN+(6.447210311379962D-03*(MIELAN^25))
2170 MIELAN=MIELAN+(5.74003766956053D-03*(MIELAN^27))
2180 MIELAN=MIELAN+(.00515330968258#*(MIELAN^29))
2190 MIELAN=MIELAN/2
2200 REM Angle of elevation is MIELAN; used formula on p. 12 of notes.
2210 REM Set shell's minimum possible height as zero feet; set gun's maximum angle of elevation as 90 degrees of arc.
2220 REM End flowchart box E3
2230 REM Begin flowchart box E4
2240 REM Set first trial angle of elevation of gun as (45 degrees of arc minus minimum angle of elevation of gun
2250 REM to reach HOSEP)/2; the reason I use 45 and not 90 degrees here is that as soon as you use an angle greater
2260 REM than 45 degrees you meet up with the possibility of the shell falling short of HOSEP.
2270 ELANG=(.785398163#-MIELAN)/2
2280 REM End flowchart box E4
2290 REM Begin flowchart box E5
2300 REM Calculate trial vertical (Z) separation resulting from ELANG and inputted values from flowchart box E2
2310 REM (except VESEP) using equation 211a; calculate difference between trial vertical separation and VESEP.
2320 A1=MV*SIN(ELANG):A2=MV*COS(ELANG)
2330 TRIZ=(HOSEP*A1)/A2
2340 TRIZ=TRIZ-G1*((HOSEP*HOSEP)/(A2*A2))
2350 REM PRINT "Trial Z and trial gun's angle of elevation are ";TRIZ,ELANG*57.29577951#:REM For debugging only
2360 REM End flowchart box E5
2370 REM Begin flowchart box E6
2380 IF ABS(TRIZ-VESEP)<.01 THEN GOTO 2470:REM That is, go to flowchart box E8. Note that this represents a 1/100
2390 REM of a foot accuracy of the shell's striking point in relation to the program's perception of the target,
2400 REM where the target is an aircraft! Compare this with the size of any conventional aircraft.
2410 REM I have tried to get the program to achieve a higher level of precision, but this only caused
2420 REM endless looping due to the problem of attainable accuracy of TRIZ and ELANG, loss of digits, etc.
2430 REM End flowchart box E6
2440 REM Begin flowchart box E7
2450 ELANG=ELANG-(((TRIZ-VESEP)/2)/CONST):GOTO 2290:REM That is, go to flowchart box E5
2460 REM End flowchart box E7
2470 REM Begin flowchart box E8
2480 REM PRINT:PRINT "Gun's angle of elevation from horizontal (depression below horizontal negative) = ";
2490 REM PRINT ELANG*57.29577951#;" decimal degrees of arc":PRINT
2500 REM Previous two lines are for debugging only.
2510 REM This next section is not in the original flowchart box E8 or any other flowchart box; however, let
2520 REM us consider this section, for all intents and purposes, to be in flowchart box E8 since its result
2530 REM (i.e. time of flight) is dependent on angle of elevation and horizontal or vertical separation (these
2540 REM last depending on the gun's angle of elevation) via a direct algebraic formula.
2550 REM Calculate the shell's time of flight.
2560 IF ELANG>1.483529864# THEN GOTO 2620:REM That is, we must use a 'special technique' to calculate
2570 REM the shell's time of flight if the gun's angle of elevation is extremely high, i.e. between
2580 REM 85 and 90 degrees of arc.
2590 TOF2=HOSEP/A2:GOTO 2680:REM Time of flight for relatively low gun angles of elevation, that is, angles of
2600 REM elevation lower than 85 degrees of arc; the GOTO etc. command in the previous line of code gets
2610 REM the program to skip over the high-gun-angle-time-of-flight algorithm.
2620 REM Calculate time of flight for gun angles of elevation between 85 and 90 degrees of arc; this uses
2630 REM equation 254a from page 254 of notes.
2640 B1=A1/G:B2=(A1*A1)/(G*G):B3=(2*VESEP)/G:TOF2=B1-(SQR(B2-B3)):REM TOF2 is time of flight in seconds
2650 REM of time for high gun angles of elevation, that is, greater than 85 degrees of arc.
2660 REM End flowchart box E8
2670 REM Begin flowchart box E9 (begin/end marker)
2680 RETURN
2690 REM End flowchart box E9
Program SPACBAT.BAS
10 REM Begin flowchart box 1
20 CLEAR:CLS:PRINT "Space battleship fire control program"
30 PRINT "for IBM PC or similar computers"
40 PRINT "which do trigonometric calculations in radians"
50 PRINT "by Ian Douglas, year 2010"
60 PRINT "This program works for any"
70 PRINT "azimuthal bearing of target ship"
80 PRINT "both at beginning and end of"
90 PRINT "shell's time of flight."
100 PRINT "Also works for any"
110 PRINT "bearing of point to aim at,"
120 PRINT "i.e. point 'C.'"
130 PRINT "X Y Z coordinates of own ship at start"
140 PRINT "of problem are 50 nautical miles by 50 nautical miles by 0."
150 PRINT "Own ship's true direction always assumed to"
160 PRINT "be a bearing of 0 degrees or 360 degrees."
170 REM This program assumes that it is working in deep space, where the shell flies straight at the target
180 REM under no gravitational influences.
190 PRINT "This program works for targets whose starting and/or ending at elevations are above or below own ship;"
200 PRINT "likewise for the shell's point of impact. If gun must be depressed to hit target, then gun's"
210 PRINT "angle of elevation is negative."
220 REM End flowchart box 1
230 REM Begin flowchart box 2
240 INPUT "Muzzle velocity in feet/sec ";MV
250 INPUT "Target spaceship's present straight line-of-sight range in feet ";TSR
260 INPUT "Target spaceship's true 'ground track' speed in feet/sec ";TSS
270 PRINT "Target spaceship's modified true plot 'ground track' direction of travel as an"
280 INPUT "azimuthal bearing in decimal degrees of arc "; TSD:TSD=TSD/57.29577951#
290 PRINT "Target spaceship's first observed angle of elevation above horizontal"
300 INPUT "in decimal degrees of arc ";TSAE:TSAE=TSAE/57.29577951#
310 PRINT "Target spaceship's rate of climb in feet/sec"
320 INPUT "(descent negative) ";TSRC
330 INPUT "Own ship's true speed in feet/sec";OSS
340 PRINT "Target spaceship's bearing off own ship"
350 INPUT "at time of gun's firing in decimal degrees of arc ";TSB:TSB=TSB/57.29577951#
360 REM End flowchart box 2
370 REM Begin flowchart box 3
380 REM Calculate target aircraft's flat range (straight line of sight range reduced to XY plane) and then
390 REM calculate beginning (i.e. initially sighted) coordinates in (feet, feet)
400 FLRANG=TSR*COS(TSAE)
410 TSPX=(FLRANG*SIN(TSB))+303805!
420 TSPY=(FLRANG*COS(TSB))+303805!
430 TSPZ=TSR*SIN(TSAE)
440 REM Above addition of 303805 is number of feet in 50 nautical mile;
450 REM recall that own ship's position is (50 naut. mile, 50 naut. mile, 0)
460 PRINT "Target spaceship's starting (i.e. initially sighted) XYZ position in (feet, feet, feet): "
470 PRINT TSPX,TSPY,TSPZ
480 REM Previous 3 lines are for debugging and error isolation only; delete or rem them out later
490 TOF1=TSR/MV
500 REM TOF1 is first approximation of shell's time of flight
510 PRINT "First approximation of shell's time of flight in seconds:"
520 PRINT TOF1
530 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
540 REM End flowchart box 3
550 REM Begin flowchart box 4
560 BEEP:REM Gets computer to beep to show that the program is looping back properly
570 REM on first (presumably wrong) guess of time of flight; need to use
580 REM BEEP since I cannot get any of my computers to 'listen' to the
590 REM scroll lock key on the keyboard(s). BEEP statement also delays operation
600 REM of program so you can see on screen whether proper reiterative looping
610 REM to get progressively better answers is taking place.
620 REM First step in this flowchart box is to get distance target spaceship
630 REM travelled during TOF (time of flight) calculated above.
640 DTST=TSS*TOF1:REM DTST is 'ground track' distance target spaceship travels during TOF1
650 REM (time of flight) calculated above
660 REM Now calculate how far own ship travelled during time of shell's flight
670 DOST=OSS*TOF1:REM DOST is 'ground track'distance own ship travelled during shell's
680 REM time of flight
690 REM Next step is to use section B of coordinate inverse program pp. 75-76
700 REM to calculate the end coordinates of the target spaceship at end of time
710 REM of flight, using azimuth of target spaceship, distance travelled, and
720 REM target spaceship's starting coordinates.
730 REM Input is target spaceship's initial sighted coordinates (TSPX, TSPY, TSPZ),
740 REM target spaceship's modified true plot direction of travel azimuthal bearing
750 REM TSD and distance target spaceship travels during shell's time of flight DTST.
760 REM Output is target spaceship's endpoint coordinates (TSEPX, TSEPY, TSEPZ).
770 REM Calculate end coordinates
780 REM of vector, from angle, distance
790 REM and starting coordinates
800 TSEPX=TSPX+(DTST*SIN(TSD)):TSEPY=TSPY+(DTST*COS(TSD)):TSEPZ=TSPZ+(TSRC*TOF1)
810 REM End calculation of target spaceship's end coordinates
820 REM Begin calculation of coordinates of point "C" (PCX, PCY, PCZ)
830 REM X coordinate of point "C" is simply x coordinate of target spaceship's end point
840 REM at end of shell's time of flight, i.e. PCX=TSEPX
850 PCX=TSEPX
860 REM Y coordinate of point "C" is the y coordinate of target spaceship's endpoint at end of time of flight of shell TSEPY
870 REM minus distance own ship travels during shell's time of flight
880 PCY=TSEPY-DOST
890 REM Z coordinate of point "C" is simply target spaceship's altitude at end of shell's time of flight, i.e. TSEPZ
900 PCZ=TSEPZ
910 REM End calculation of coordinates of point "C"
920 PRINT "'Ground track' distance target spaceship travels during shell's time of flight (feet):"
930 PRINT DTST
940 PRINT "'Ground track' distance own ship travels during shell's time of flight (feet):"
950 PRINT DOST
960 PRINT "Target spaceship's XYZ end coordinates (feet, feet, feet):"
970 PRINT TSEPX,TSEPY, TSEPZ
980 PRINT "XYZ coordinates of point 'C' (feet, feet, feet):"
990 PRINT PCX,PCY,PCZ
1000 REM Previous 8 lines are for debugging and error isolation only; delete or rem them out later
1010 REM End flowchart box 4
1020 REM Begin flowchart box 5
1030 REM Calculate straight line distance between own ship at time of gun's firing
1040 REM and point "C," DOSPC
1050 REM Recall that own ship's coordinates at time of gun's firing are
1060 REM (303805 feet, 303805 feet, 0 feet) (i.e. 50 nautical miles by 50 nautical miles by 0)
1070 DOSPC=SQR(((PCX-303805!)^2)+((PCY-303805!)^2)+(PCZ^2))
1080 PRINT "Distance between own ship and point 'C' when gun is fired (feet):"
1090 PRINT DOSPC
1100 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
1110 REM End flowchart box 5
1120 REM Begin flowchart box 6
1130 TOF2=DOSPC/MV
1140 PRINT "New time of flight (seconds of time):"
1150 PRINT TOF2
1160 REM Previous 2 lines are for debugging and error isolation only; delete or rem them out later
1170 REM End flowchart box 6
1180 REM Begin flowchart box 7
1190 IF (ABS(TOF2-TOF1))<.001 THEN GOTO 1280
1200 REM Note that the above line of code gives 1/1000 of a second accuracy!! I have found through
1210 REM experimentation that seeking higher accuracy will in some cases cause
1220 REM infinite looping of the program in a search for unnecessarily high accuracy.
1230 REM It is debatable whether in this program you need even 1/100 of a second accuracy.
1240 REM End flowchart box 7
1250 REM Begin flowchart box 8
1260 TOF1=TOF2:GOTO 550
1270 REM End flowchart box 8
1280 REM Begin flowchart box 9
1290 REM End flowchart box 9
1300 REM Begin flowchart box 10
1310 REM Calculate gun's bearing (not angle of elevation yet) using only X Y coordinates of point "C"
1320 REM From flowchart box 4, and using coordinates of own ship at time of gun's firing;
1330 REM Convert points to vector
1340 REM 1st point is coordinates of own ship at gun's time of firing (303805 feet, 303805 feet)
1350 REM 2nd point is coordinates of point "C," (PCX, PCY)
1360 DXAA=PCX-303805!:DYAA=PCY-303805!
1370 IF DXAA=0 THEN GOTO 1430
1380 REM Above line prevents division by zero in next line
1390 SLOPE=DYAA/DXAA:REM Slope
1400 ANGLE#=ATN(SLOPE):REM Angle referred to standard position, calculated from slope
1410 IF DXAA<0 THEN ANGLE#=ANGLE#+3.141592654#
1420 IF DXAA>0 AND DYAA<0 THEN ANGLE#=ANGLE#+6.283185307#
1430 IF DXAA=0 AND DYAA=0 THEN ANGLE#=0
1440 IF DXAA=0 AND DYAA>0 THEN ANGLE#=1.570796327#
1450 IF DXAA<0 AND DYAA=0 THEN ANGLE#=3.141592654#
1460 IF DXAA=0 AND DYAA<0 THEN ANGLE#=4.712388981#
1470 IF DXAA>0 AND DYAA=0 THEN ANGLE#=6.283185307#
1480 GUNAZ#=6.283185307#-ANGLE#+1.570796327#:REM 360 degrees minus angle plus 90 degrees
1490 IF GUNAZ#>6.283185307# THEN GUNAZ#=GUNAZ#-6.283185307#
1500 REM Calculate distance between own ship and 'ground track' position of point 'C,' that is,
1510 REM 'top view' position of point 'C'
1520 DOCG=SQR((DXAA^2)+(DYAA^2)):REM Above distance is calculated via Pythagoras
1530 REM Calculate gun's angle of elevation using arctan of (PCZ/DOCG)
1540 THETA1=ATN(PCZ/DOCG)
1550 REM THETA1 is the gun's angle of elevation such that the gun points directly at point "C"
1560 REM when it is fired
1570 REM End flowchart box 10
1580 REM Begin flowchart box 11
1590 REM Gun's angle of elevation is previously calculated THETA1
1600 PRINT "Gun's azimuth to shoot along: ";GUNAZ#*57.29577951#;" decimal degrees of arc"
1610 PRINT "Gun's angle of elevation: "; THETA1*57.29577951#;" decimal degrees of arc"
1620 REM If latest value for time of flight is to be printed here, that value is TOF2.
1630 IF TOF2<301 THEN GOTO 1680
1640 PRINT:PRINT "Time of shell's flight is greater than 5 minutes; target spaceship"
1650 PRINT "is effectively out of range.":PRINT
1660 REM End flowchart box 11
1670 REM Begin flowchart box 12
1680 END
1690 REM End flowchart box 12
Appendix, January 2012: Is an aircraft or airborne target within range of an antiaircraft gun?
A surface-to-surface artillery piece is surrounded by a 2 dimensional "circle of death" whose radius is equivalent
to the gun's surface (of the ground or sea) range; this circle's plane is also on the ground or sea surface, level with the gun. Analysis using basic trajectory equations will show that the space within range of an antiaircraft gun's upward-only
and upward-to-peak trajectories is an "ellipsoid of death," that is, a (3 dimensional) bubble-like ellipsoid of
revolution. The semimajor axis of this ellipsoid is horizontal, and equivalent to half the gun's ground/sea range.
The minor axis is vertical, and is equivalent to the shell's maximum vertically-fired attainable height. The gun's
location just touches the bottom of the ellipsoid, at the bottom of the minor axis; the major axis's closest point
to the gun is above the gun by a distance equal to half the shell's maximum attainable vertically-fired height.
The major axis is exactly twice as long as the minor axis, and the "bubble" space of the ellipsoid of death has a
flattened appearance. All the above statements hold true given assumptions I have made elsewhere in this article,
namely, that there is no air resistance, that the earth is flat, gravity is the only force other than the gun's
initial impulse in effect, the earth is not rotating, that there is no wind (or even air for that matter), that
there are no Magnus or Coriolis effects, etc.
I have added some program listings to this article to determine whether a target is within range of an antiaircraft gun. The geometry of ellipses (related to these types of problems) is discussed (besides many other sources) in Krafft Ehricke's book "Space Flight Vol. I, Environment and Celestial Mechanics" (Princeton, N.J: D. Van Nostrand Company, Inc., 1960, pp. 293 et seq.). At a future date, I hope to integrate such a program with UPART.BAS and AIRBAT.BAS. As the reader may have noticed by now, my programs UPART.BAS and AIRBAT.BAS do not indicate or trap situations where the target aircraft is out of range; in my programs at present, these situations may cause infinite looping or an alarm-like exit to the operating system.
Below (program IOELLIP.BAS) is a program to determine whether a 2 dimensional point is inside or outside of a 2D ellipse whose major axis is on the x axis. Part of the algorithm used here is from Ehricke's book mentioned above. This program would be useful for problems where the target aircraft is flying on a straight ground track, this ground track passing directly over the site of the gun. Regarding coordinate systems, there are two in use here: the gun's coordinate system (implied in this situation but not used in the program below) which has the location of the gun at the origin; and the ellipse's coordinate system, which has the centre of the ellipse at its origin. In a combat situation, the gun and ellipse centre are separated vertically; the ellipse centre and origin are above the gun location and origin by a distance equal to the ellipse's semi minor axis length. That semi minor axis length is in turn equal to half the shell's maximum attainable vertically-fired height. So ultimately we would have to convert between coordinate systems thus:
YE = ellipse Y coordinate
YG = gun Y coordinate
B = ellipse semi minor axis length (i.e. half the shell's maximum attainable vertically-fired height)
A gun Y coordinate is always positive in these types of problem.
To convert from YE to YG:
if YE is positive, then YG = ABS(B) + YE
if YE is negative, then YG = ABS(B) - ABS(YE) (works only when gun Y coordinate is positive, which applies to our problem)
[ ABS (expression/variable/number) means "the absolute value of (expression/variable/number) ]
To convert from YG to YE:
If YG < ABS(B) then YE = - {[ABS(B)] - YG}
If YG > ABS(B) then YE = YG - ABS(B)
Program IOELLIP.BAS
10 REM Begin flowchart box 1
20 CLS:CLEAR:PRINT "Determination of whether a given xy point in the 2D"
30 PRINT "Cartesian plane is inside or outside an ellipse"
40 PRINT "of given major and minor axis lengths"
50 PRINT "with the ellipse centre being at (0,0)."
60 PRINT
70 REM End flowchart box 1
80 REM Begin flowchart box 2
90 PRINT "Enter the data for the ellipse."
100 INPUT "Semi major axis length ";MAJA
110 INPUT "Semi minor axis length ";MINA
120 PRINT:PRINT "Enter the coordinates for the point whose general location (i.e. inside or outside"
130 PRINT "the ellipse) we wish to determine."
140 PRINT:INPUT "Point's x coordinate ";XB
150 INPUT "Point's y coordinate ";YB
160 REM End flowchart box 2
170 REM Begin flowchart box 3
180 XA=0:YA=0:REM Note that the starting coordinates XA and YA are at the centre of the ellipse
190 REM which is also the origin (0,0).
200 DX=XB-XA:DY=YB-YA
210 D=SQR((DX^2)+(DY^2)):REM Distance
220 IF DX=0 THEN GOTO 280
230 REM Above line prevents division by zero in next line
240 S=DY/DX:REM Slope
250 AN=ATN(S)
260 IF DX<0 THEN AN=AN+3.141592654#
270 IF DX>0 AND DY<0 THEN AN=AN+6.283185307#
280 IF DX=0 AND DY=0 THEN AN=0
290 IF DX=0 AND DY>0 THEN AN=1.570796327#
300 IF DX<0 AND DY=0 THEN AN=3.141592654#
310 IF DX=0 AND DY<0 THEN AN=4.712388981#
320 IF DX>0 AND DY=0 THEN AN=6.283185307#
330 REM End flowchart box 3
340 REM Begin flowchart box 4
350 REM Use Ehricke Vol. 1 p. 293 eq. 4-92.
360 REM Plug in central polar anomaly 'psi'
370 REM (which is the variable AN obtained in flowchart
380 REM box 3) to get the central polar radius vector
390 REM of the ellipse 'rho.' Above mentioned
400 REM central polar anomaly is the same for both
410 REM the ellipse and the given point.
420 A2=MAJA^2:B2=MINA^2:COS2P=(COS(AN))^2:TERM1=A2*B2:TERM2=(A2-B2)*COS2P
430 TERM3=A2-TERM2:REM Must keep this line and the next two non-remmed-out lines separate from one another
440 REM or program will not work! Logically they should work but the interpreter will not recognize this!
450 IF TERM3=0 THEN TERM3=.0000001
460 RHO=TERM1/TERM3:RHO=SQR(RHO)
470 REM Line before previous line prevents division by zero later in that same line to obtain 'rho'.
480 PRINT "Rho (i.e. central polar radius vector of"
490 PRINT "the ellipse for the same central"
500 PRINT "polar anomaly as the given point) ";RHO
510 PRINT "Distance of given point from centre of"
520 PRINT "ellipse ";D
530 PRINT "Psi (i.e. central polar anomaly of the point and of the ellipse intersection point of that line) "
540 PRINT AN;" radians":PRINT
550 REM End flowchart box 4
560 REM Begin flowchart box 5 (off page connector)
570 REM
580 REM End flowchart box 5 (off page connector)
590 REM Begin flowchart box 6 (off page connector)
600 REM
610 REM End flowchart box 6 (off page connector)
620 REM Begin flowchart box 7
630 IF RHO>D THEN PRINT "Given point is inside ellipse ":GOTO 660:REM Latter part of this line is flowchart box 8
640 PRINT "Given point is on or outside ellipse":REM This is flowchart box 9
650 REM End flowchart box 7
660 REM Begin flowchart box 10
670 END
680 REM End flowchart box 10
Here is a program to determine whether a 3 dimensional XYZ point is inside or outside an ellipsoid of revolution, sort of a "geodesic ellipsoid." The program below works with the left hand orthogonal coordinate system, but probably works equally well with the right hand orthogonal system. The coordinate system being used must have Z being vertical, with the semi minor axis being on Z. The ellipsoid is like two frisbees glued together by their bottom surfaces, clamshell-like or flying-saucer-like; its major axes are on X and Y, and are both the same length. The ellipsoid used is not a "rugby football" or "torpedo" ellipsoid, with differing lengths for each ellipsoid axis.
Program IOELOID.BAS
10 CLS:CLEAR:PRINT "Program to determine whether a given XYZ point"
20 PRINT "is inside or outside an ellipsoid with semi minor axis being on the Z axis"
30 PRINT "and semi major axes, both of same length, being on X and Y axes."
40 PRINT "Note that this is a 'clamshell flying saucer' ellipsoid and not a 'torpedo' or 'rugby football' ellipsoid."
50 PRINT "Program calculates central polar latitude and radius vector"
60 PRINT "given a point's XYZ coordinates and lengths of"
70 PRINT "major and minor axes. Coordinate system"
80 PRINT "used here is right handed, sort of a 'geodesic ellipsoid' with"
90 PRINT "'north pole' being on positive Z axis"
100 PRINT "and X and Y being 'equatorial radii.'"
110 REM For in-or-out-of-flat-ellipse calculations, see Ehricke, "Space Flight Vol. 1", page 293, eq. 4-92
120 PRINT:PRINT "Enter the given point's XYZ coordinates"
130 PRINT:INPUT "X = ";GPX
140 INPUT "Y = ";GPY
150 INPUT "Z = ";GPZ
160 PRINT:PRINT "Enter the ellipsoid's semi major and semi minor axis lengths."
170 PRINT:INPUT "Semi major axis length = ";SMAJA
180 INPUT "Semi minor axis length = ";SMINA
190 REM Calculate the non-reduced non-projected horizontal distance from the origin to the
200 REM given point (i.e. to the location of the given point as seen projected onto the XY plane as seen from
210 REM well above the ellipsoid's "north pole" (the "north pole" being the Z axis).
220 HDIST=SQR((GPX^2)+(GPY^2))
230 REM Calculate the straight line distance from
240 REM the origin to the given point
250 SLD=SQR((GPX^2)+(GPY^2)+(GPZ^2))
260 REM Calculate the central polar non-reduced non-projected anomaly "psi" from the XY plane
270 REM to the given point; this is not a full circle arctan situation, as the denominator argument
280 REM is always positive.
290 DENOM1=SQR((GPX^2)+(GPY^2)):IF DENOM1=0 THEN DENOM1=.0000001:REM Last statement prevents division by zero in next line.
300 PSI=ATN(GPZ/DENOM1)
310 REM Calculate central polar radius vector of ellipsoid
320 TERM1=SMAJA^2:TERM2=SMINA^2:TERM3=TERM1*TERM2:TERM4=TERM1-TERM2
330 TERM5=(COS(PSI))^2:TERM6=TERM4*TERM5
340 TERM7=TERM1-TERM6:REM Must keep this line and the next two non-remmed-out lines separate from one another
350 REM or program will not work! Logically they should work but the interpreter will not recognize this!
360 IF TERM7=0 THEN TERM7=.0000001
370 TERM8=TERM3/TERM7:RHO=SQR(TERM8)
380 REM Line before previous line prevents division by zero later in that same line.
390 PRINT:PRINT "Straight line distance from origin to given point: "
400 PRINT SLD
410 PRINT:PRINT "Central polar radius vector from origin to ellipsoid:"
420 PRINT RHO:PRINT
430 IF RHO>SLD THEN PRINT "Given point is inside ellipsoid":GOTO 450
440 PRINT "Given point is on or outside ellipsoid"
450 PRINT:END
The following two programs (GUNEOD1.BAS and EIDGUN2.BAS) deal with the problem of taking a stationary airborne
target's antiaircraft-gun-centred coordinates and converting them to ellipsoid-of-death-centred coordinates; then
the programs determine whether the airborne target is within range of the antiaircraft gun. Upward-only trajectories
are considered, not upward-to-peak-to-downward trajectories. As before, the only forces assumed to be acting on the
shell are the gun's initial impulse and gravity. Gun-centred coordinates have the gun's location at the origin; ellipse
of death centred coordinates, or ellipsoid of death centred coordinates have the centre of the ellipse or ellipsoid
at the origin. GUNEOD1.BAS solves the problem in 2 dimensions, where x is horizontal and y is vertical; EIDGUN2.BAS
solves the problem in 3 dimensions, where x and y are horizontal and z is vertical. In the 2D version of the
problem, x is the same for both gun-centred and ellipse of death centred coordinates; the only conversion that
needs to be made is in y (vertical). In the 3D version of the problem, x and y are the same for both gun-centred
and ellipse of death centred coordinates; the only conversion that needs to be made is in z (vertical).
Program GUNEOD1.BAS
10 REM Begin flowchart box 1
20 CLS:CLEAR:PRINT "Determination of whether a target's xy location in the 2D"
30 PRINT "Cartesian plane is inside or outside an AA gun's ellipse of death (EOD);"
40 PRINT "the 2D plane has its y axis as vertical. Gun-centred coordinates used below have the location"
50 PRINT "of the gun at the origin; EOD-centred coordinates used below have the centre of the ellipse at"
60 PRINT "the origin."
70 PRINT
80 REM End flowchart box 1
90 REM Begin flowchart box 2
100 PRINT "Enter the data for the ellipse."
110 INPUT "Semi major axis length ";MAJA
120 INPUT "Semi minor axis length ";MINA
130 PRINT:PRINT "Enter the coordinates for the target's point whose general location (i.e. inside or outside"
140 PRINT "the ellipse) we wish to determine."
150 PRINT:INPUT "Target's gun-centred x coordinate ";XBG:XB=XBG:REM Last statement is simply intended to keep program's
160 REM and logic easier to follow than would be otherwise.
170 INPUT "Target's gun-centred y coordinate ";YBG
180 REM Convert y (vertical) coordinate from gun centred coordinate to ellipse-of-death-centred coordinate
190 IF YBG=ABS(MINA) THEN YB=YBG-ABS(MINA)
210 IF YBG=0 THEN YB=-(ABS(MINA))
220 REM XB, YB obtained above are EOD-centred x, y coordinates.
230 PRINT "Target's EOD-centred x coordinate ";XB
240 PRINT "Target's EOD-centred y coordinate ";YB
250 REM End flowchart box 2
260 REM Begin flowchart box 3
270 XA=0:YA=0:REM Note that the starting coordinates XA and YA are at the centre of the ellipse
280 REM which is also the origin (0,0).
290 DX=XB-XA:DY=YB-YA
300 D=SQR((DX^2)+(DY^2)):REM Distance
310 IF DX=0 THEN GOTO 370
320 REM Above line prevents division by zero in next line
330 S=DY/DX:REM Slope
340 AN=ATN(S)
350 IF DX<0 THEN AN=AN+3.141592654#
360 IF DX>0 AND DY<0 THEN AN=AN+6.283185307#
370 IF DX=0 AND DY=0 THEN AN=0
380 IF DX=0 AND DY>0 THEN AN=1.570796327#
390 IF DX<0 AND DY=0 THEN AN=3.141592654#
400 IF DX=0 AND DY<0 THEN AN=4.712388981#
410 IF DX>0 AND DY=0 THEN AN=6.283185307#
420 REM End flowchart box 3
430 REM Begin flowchart box 4
440 REM Use Ehricke Vol. 1 p. 293 eq. 4-92.
450 REM Plug in central polar anomaly 'psi'
460 REM (which is the variable AN obtained in flowchart
470 REM box 3) to get the central polar radius vector
480 REM of the ellipse 'rho.' Above mentioned
490 REM central polar anomaly is the same for both
500 REM the ellipse and the target.
510 A2=MAJA^2:B2=MINA^2:COS2P=(COS(AN))^2:TERM1=A2*B2:TERM2=(A2-B2)*COS2P
520 TERM3=A2-TERM2:REM Must keep this line and the next two non-remmed-out lines separate from one another
530 REM or program will not work! Logically they should work but the interpreter will not recognize this!
540 IF TERM3=0 THEN TERM3=.0000001
550 RHO=TERM1/TERM3:RHO=SQR(RHO)
560 REM Line before previous line prevents division by zero later in that same line to obtain 'rho'.
570 PRINT "Rho (i.e. central polar radius vector of"
580 PRINT "the ellipse for the same central"
590 PRINT "polar anomaly as the target) ";RHO
600 PRINT "Distance of target from centre of"
610 PRINT "ellipse ";D
620 PRINT "Psi (i.e. central polar anomaly of the point and of the ellipse intersection point of that line) "
630 PRINT AN;" radians":PRINT
640 REM End flowchart box 4
650 REM Begin flowchart box 5 (off page connector)
660 REM
670 REM End flowchart box 5 (off page connector)
680 REM Begin flowchart box 6 (off page connector)
690 REM
700 REM End flowchart box 6 (off page connector)
710 REM Begin flowchart box 7
720 IF RHO>D THEN PRINT "Target is inside ellipse of death":GOTO 750:REM Latter part of this line is flowchart box 8
730 PRINT "Target is on or outside ellipse of death":REM This is flowchart box 9
740 REM End flowchart box 7
750 REM Begin flowchart box 10
760 END
770 REM End flowchart box 10
Program EIDGUN2.BAS
10 CLS:CLEAR:PRINT "Program to determine whether a target's XYZ is inside or outside an"
20 PRINT "AA gun's ellipsoid of death (EOD) with semi minor axis being on the Z axis"
30 PRINT "and semi major axes, both of same length, being on X and Y axes."
40 PRINT "Note that this is a 'clamshell flying saucer' ellipsoid and not a 'torpedo' or 'rugby football' ellipsoid."
50 PRINT "Program calculates central polar latitude and radius vector"
60 PRINT "given a target's XYZ coordinates and lengths of major and minor axes. Coordinate system"
70 PRINT "used here is right handed, sort of a 'geodesic ellipsoid' with 'north pole' being on positive Z axis"
80 PRINT "and X and Y being 'equatorial radii.' Gun-centred coordinates used below have the location"
90 PRINT "of the gun at the origin; EOD-centred coordinates used below have the centre of the EOD"
100 PRINT "at the origin."
110 REM For in-or-out-of-flat-ellipse calculations, see Ehricke, "Space Flight Vol. 1", page 293, eq. 4-92
120 PRINT:PRINT "Enter the target's gun-centred XYZ coordinates"
130 PRINT:INPUT "X = ";GPGCX
140 INPUT "Y = ";GPGCY
150 INPUT "Z = ";GPGCZ
160 GPX=GPGCX:GPY=GPGCY:REM Last statement is simply intended to keep program's logic easier to follow
170 REM than would be otherwise.
180 PRINT:PRINT "Enter the ellipsoid's semi major and semi minor axis lengths."
190 PRINT:INPUT "Semi major axis length = ";SMAJA
200 INPUT "Semi minor axis length = ";SMINA
210 REM Convert z (vertical) coordinate from gun centred coordinate to EOD-centred coordinate
220 IF GPGCZ=ABS(SMINA) THEN GPZ=GPGCZ-ABS(SMINA)
240 IF GPGCZ=0 THEN GPZ=-(ABS(SMINA))
250 REM GPX, GPY, GPZ obtained above are EOD-centred x, y, z coordinates.
260 PRINT:PRINT "Target's EOD-centred x, y, z coordinates:"
270 PRINT:PRINT "X = ";GPX
280 PRINT "Y = ";GPY
290 PRINT "Z = ";GPZ
300 REM Calculate the non-reduced non-projected horizontal distance from the origin to the
310 REM target (i.e. to the location of the target as seen projected onto the XY plane as seen from
320 REM well above the ellipsoid's "north pole") (the "north pole" being the Z axis).
330 HDIST=SQR((GPX^2)+(GPY^2))
340 REM Calculate the straight line distance from
350 REM the origin to the target
360 SLD=SQR((GPX^2)+(GPY^2)+(GPZ^2))
370 REM Calculate the central polar non-reduced non-projected anomaly "psi" from the XY plane
380 REM to the target; this is not a full circle arctan situation, as the denominator argument
390 REM is always zero or positive.
400 DENOM1=SQR((GPX^2)+(GPY^2))
410 IF DENOM1=0 THEN DENOM1=.0000001:REM Last statement prevents division by zero in next line.
420 PSI=ATN(GPZ/DENOM1)
430 REM Calculate central polar radius vector of ellipsoid
440 TERM1=SMAJA^2:TERM2=SMINA^2:TERM3=TERM1*TERM2:TERM4=TERM1-TERM2
450 TERM5=(COS(PSI))^2:TERM6=TERM4*TERM5
460 TERM7=TERM1-TERM6:REM Must keep this line and the next two non-remmed-out lines separate from one another
470 REM or program will not work! Logically they should work but the interpreter will not recognize this!
480 IF TERM7=0 THEN TERM7=.0000001
490 TERM8=TERM3/TERM7:RHO=SQR(TERM8)
500 REM Line before previous line prevents division by zero later in that same line.
510 PRINT:PRINT "Straight line distance from origin to target: "
520 PRINT SLD
530 PRINT:PRINT "Central polar radius vector from origin to ellipsoid:"
540 PRINT RHO:PRINT
550 IF RHO>SLD THEN PRINT "Target is inside ellipsoid of death":GOTO 570
560 PRINT "Target is on or outside ellipsoid of death"
570 PRINT:END
Program UPARTR.BAS :This program expands on program UPART.BAS so as to trap instances where the target is outside
the EOD and therefore out of range. Also note that this program has a weakness also present in UPART.BAS, which is
that it goes into an infinite loop if the fire and hit solution involves the gun having an angle of elevation of
greater than about 78 degrees: see the introductory notes for the program UPART.BAS above for a fuller explanation.
This program appears to solve the problem in UPART.BAS of the program simply crashing, or looping forever, when
the target is out of range.
10 REM Begin flowchart box E1 (begin/end marker)
20 REM End flowchart box E1; the flowchart numbers in this program correspond to flowchart on p. 251 of notes.
30 REM Begin flowchart box E2
40 PRINT "Program to use a variant of successive approximation to get angle of elevation"
50 PRINT "of a gun firing at a target separated from the gun significantly"
60 PRINT "separated from it both vertically and horizontally (e.g. an aircraft) using equation"
70 PRINT "211a on page 211 of notes; by Ian Douglas, year 2010. This program also works"
80 PRINT "for targets below gun. This program is for upward-only trajectories if target is above gun"
90 PRINT "and not for upward-to-peak-to-downward trajectories."
100 REM Program is intended to make the problem manageable by taking the maximum possible shell height the gun
110 REM can attain and obtaining a constant to multiply 90 degrees of arc (1.570796327 radians) (gun's maximum angle of
120 REM elevation) by to get any given possible height the shell can attain; this enables a resultant vertical
130 REM separation of the target from gun to be compared to an actual vertical separation attained by a trial angle
140 REM of elevation of gun, for comparison purposes; also this enables the successive approximation process to occur.
150 REM For this reason, I do not recommend using this program to compute data concerning guns with greater
160 REM than 10,000 feet per second muzzle velocities; but even the most powerful guns have about 1/4 that.
170 REM Both the program and its subroutine GUNEOD2 assume that the problem is in a 2 dimensional plane, with
180 REM the plane's horizontal portion being on the line on the (perfectly flat and level) ground separating
190 REM the target from the gun horizontally, and the vertical portion of the plane being the vertical separation
200 REM of the gun from the target.
210 REM Also note that the program will go into an infinite loop if the target is directly above gun, even though
220 REM a fire and hit solution is possible in such a case.
230 PRINT "Vertical separation (Y) of target from gun (feet) "
240 PRINT "(negative for a target below gun)"
250 INPUT VESEP
260 INPUT "Horizontal separation (in XY plane) of target from gun (feet) ";HOSEP
270 INPUT "Gun's muzzle velocity (feet per second) ";MV
280 G=32:G1=G/2:REM Acceleration due to gravity in feet per second per second; also half this for use in
290 REM instances of equation 211a.
300 REM End flowchart box E2
310 REM Begin flowchart box E3
320 REM Calculate the shell's maximum attainable height.
330 MAXHT=(MV*MV)/(2*G)
340 REM Calculate the semi major axis length of the Ellipse of Death (EOD); this is half the gun's surface-to-
350 REM surface range.
360 MAJA=((MV*MV)/G)/2
370 REM Calculate the semi minor axis length of the Ellipse of Death (EOD); this is half the shell's maximum
380 REM attainable height, i.e. the maximum height attained by the shell when fired straight up.
390 MINA=MAXHT/2:GOSUB 970
400 REM Calculate the constant mentioned in remarks for flowchart box E2.
410 CONST=MAXHT/1.570796327#
420 REM Calculate the minimum angle of elevation of gun needed to reach HOSEP.
430 MIELAN=(HOSEP*G)/(MV*MV)
440 REM The following is an arc sine routine with the final angle divided by 2 for the purpose of the
450 REM above formula to get the minimum necessary angle of elevation given HOSEP, muzzle velocity, and G.
460 MIELAN=MIELAN+((MIELAN^3)/6)+(.075*(MIELAN^5))
470 MIELAN=MIELAN+(4.464285714285714D-02*(MIELAN^7))
480 MIELAN=MIELAN+(3.038194444444444D-02*(MIELAN^9))
490 MIELAN=MIELAN+(2.237215909090909D-02*(MIELAN^11))
500 MIELAN=MIELAN+(1.735276442307692D-02*(MIELAN^13))
510 MIELAN=MIELAN+(.01396484375#*(MIELAN^15))
520 MIELAN=MIELAN+(1.155180089613971D-02*(MIELAN^17))
530 MIELAN=MIELAN+(9.761609529194079D-03*(MIELAN^19))
540 MIELAN=MIELAN+(8.390335809616816D-03*(MIELAN^21))
550 MIELAN=MIELAN+(7.312525877035907D-03*(MIELAN^23))
560 MIELAN=MIELAN+(6.447210311379962D-03*(MIELAN^25))
570 MIELAN=MIELAN+(5.74003766956053D-03*(MIELAN^27))
580 MIELAN=MIELAN+(.00515330968258#*(MIELAN^29))
590 MIELAN=MIELAN/2
600 REM Angle of elevation is MIELAN; used formula on p. 12 of notes.
610 REM Set shell's minimum possible height as zero feet; set gun's maximum angle of elevation as 90 degrees of arc.
620 REM End flowchart box E3
630 REM Begin flowchart box E4
640 REM Set first trial angle of elevation of gun as (45 degrees of arc minus minimum angle of elevation of gun
650 REM to reach HOSEP)/2; the reason I use 45 and not 90 degrees here is that as soon as you use an angle greater
660 REM than 45 degrees you meet up with the possibility of the shell falling short of HOSEP.
670 ELANG=(.785398163#-MIELAN)/2
680 REM End flowchart box E4
690 REM Begin flowchart box E5
700 REM Calculate trial vertical (Y) separation resulting from ELANG and inputted values from flowchart box E2
710 REM (except VESEP) using equation 211a; calculate difference between trial vertical separation and VESEP.
720 A1=MV*SIN(ELANG):A2=MV*COS(ELANG)
730 TRIZ=(HOSEP*A1)/A2
740 TRIZ=TRIZ-G1*((HOSEP*HOSEP)/(A2*A2))
750 REM PRINT "Trial Z and trial gun's angle of elevation are ";TRIZ,ELANG*57.29577951#:REM For debugging only
760 REM End flowchart box E5
770 REM Begin flowchart box E6
780 IF ABS(TRIZ-VESEP)<.01 THEN GOTO 870:REM That is, go to flowchart box E8. Note that this represents a 1/100
790 REM of a foot accuracy of the shell's striking point in relation to the program's perception of the target,
800 REM where the target is an aircraft! Compare this with the size of any conventional aircraft.
810 REM I have tried to get the program to achieve a higher level of precision, but this only caused
820 REM endless looping due to the problem of attainable accuracy of TRIZ and ELANG, loss of digits, etc.
830 REM End flowchart box E6
840 REM Begin flowchart box E7
850 ELANG=ELANG-(((TRIZ-VESEP)/2)/CONST):GOTO 690:REM That is, go to flowchart box E5
860 REM End flowchart box E7
870 REM Begin flowchart box E8
880 PRINT:PRINT "Gun's angle of elevation from horizontal (depression below horizontal negative) = ";
890 PRINT ELANG*57.29577951#;" decimal degrees of arc":PRINT
900 REM End flowchart box E8
910 REM Begin flowchart box E9 (begin/end marker)
920 END
930 REM End flowchart box E9
940 REM
950 REM * * * Subroutines * * *
960 REM
970 REM Subroutine GUNEOD2, derived from program GUNEOD1.BAS; flowchart box numbers in this subroutine
980 REM refer to notes about GUNEOD1.BAS.
990 REM Begin flowchart box 1
1000 REM Determination of whether a target's xy location in the 2D
1010 REM Cartesian plane is inside or outside an AA gun's ellipse of death (EOD);
1020 REM the 2D plane has its y axis as vertical. Gun-centred coordinates used below have the location
1030 REM of the gun at the origin; EOD-centred coordinates used below have the centre of the ellipse at
1040 REM the origin. Enter the subroutine with the EOD's semi major axis length MAJA, semi minor axis length MINA,
1050 REM horizontal separation of target from gun HOSEP, and vertical separation of target from gun VESEP. You
1060 REM leave the program knowing whether the target is within range or out of range. If the target is within
1070 REM range, control goes back to the main program. If the target is out of range, this subroutine says so
1080 REM and terminates the program.
1090 PRINT
1100 REM End flowchart box 1
1110 REM Begin flowchart box 2
1120 XBG=HOSEP:YBG=VESEP:XB=XBG:REM Last statement is simply intended to keep program's
1130 REM and logic easier to follow than would be otherwise.
1140 REM Convert y (vertical) coordinate from gun centred coordinate to ellipse-of-death-centred coordinate
1150 IF YBG=ABS(MINA) THEN YB=YBG-ABS(MINA)
1170 IF YBG=0 THEN YB=-(ABS(MINA))
1180 REM XB, YB obtained above are EOD-centred x, y coordinates.
1190 PRINT "Target's EOD-centred x coordinate ";XB
1200 PRINT "Target's EOD-centred y coordinate ";YB
1210 REM End flowchart box 2
1220 REM Begin flowchart box 3
1230 XA=0:YA=0:REM Note that the starting coordinates XA and YA are at the centre of the ellipse
1240 REM which is also the origin (0,0).
1250 DX=XB-XA:DY=YB-YA
1260 D=SQR((DX^2)+(DY^2)):REM Distance
1270 IF DX=0 THEN GOTO 1330
1280 REM Above line prevents division by zero in next line
1290 S=DY/DX:REM Slope
1300 AN=ATN(S)
1310 IF DX<0 THEN AN=AN+3.141592654#
1320 IF DX>0 AND DY<0 THEN AN=AN+6.283185307#
1330 IF DX=0 AND DY=0 THEN AN=0
1340 IF DX=0 AND DY>0 THEN AN=1.570796327#
1350 IF DX<0 AND DY=0 THEN AN=3.141592654#
1360 IF DX=0 AND DY<0 THEN AN=4.712388981#
1370 IF DX>0 AND DY=0 THEN AN=6.283185307#
1380 REM End flowchart box 3
1390 REM Begin flowchart box 4
1400 REM Use Ehricke Vol. 1 p. 293 eq. 4-92.
1410 REM Plug in central polar anomaly 'psi'
1420 REM (which is the variable AN obtained in flowchart
1430 REM box 3) to get the central polar radius vector
1440 REM of the ellipse 'rho.' Above mentioned
1450 REM central polar anomaly is the same for both
1460 REM the ellipse and the target.
1470 A2=MAJA^2:B2=MINA^2:COS2P=(COS(AN))^2:TERM1=A2*B2:TERM2=(A2-B2)*COS2P
1480 TERM3=A2-TERM2:REM Must keep this line and the next two non-remmed-out lines separate from one another
1490 REM or program will not work! Logically they should work but the interpreter will not recognize this!
1500 IF TERM3=0 THEN TERM3=.0000001
1510 RHO=TERM1/TERM3:RHO=SQR(RHO)
1520 REM Line before previous line prevents division by zero later in that same line to obtain 'rho'.
1530 PRINT "Rho (i.e. central polar radius vector of"
1540 PRINT "the ellipse for the same central"
1550 PRINT "polar anomaly as the target) ";RHO
1560 PRINT "Distance of target from centre of"
1570 PRINT "ellipse ";D
1580 PRINT "Psi (i.e. central polar anomaly of the point and of the ellipse intersection point of that line) "
1590 PRINT AN;" radians":PRINT
1600 REM End flowchart box 4
1610 REM Begin flowchart box 5 (off page connector)
1620 REM
1630 REM End flowchart box 5 (off page connector)
1640 REM Begin flowchart box 6 (off page connector)
1650 REM
1660 REM End flowchart box 6 (off page connector)
1670 REM Begin flowchart box 7
1680 PRINT "EOD semimajor axis length = ";MAJA;" EOD semiminor axis = ";MINA
1690 IF RHO>D THEN PRINT "Target is inside ellipse of death and therefore within gun's range":RETURN
1700 REM Latter part of above line is similar to flowchart box 8
1710 PRINT "Target is on or outside ellipse of death and therefore outside of gun's range"
1720 REM Above line is flowchart box 9
1730 REM End flowchart box 7
1740 REM Begin flowchart box 10
1750 END
1760 REM End flowchart box 10