10 REM PROGRAM [MCTRAJ.BAS], MAR. 1987. [TANDY 1000 HX, 07/90, REV. 02/93]
20 CLS
30 KEY OFF
40 COLOR 7, 1, 8: CLS
50 KEY ON
60 DEFDBL A-H, M, O-Z
70 REM POINT MASS TRAJECTORIES FOR SMALL ARMS
80 REM THE PROGRAM REQUIRES AN INPUT TABLE OF
90 REM DRAG COEFFICIENT (CD) VERSUS MACH NUMBER (M).
100 REM ADDITIONAL REQUIRED INPUTS ARE:
110 REM MUZZLE VELOCITY (FT/SEC); BALLISTIC COEFFICIENT (LB/IN 2);
120 REM HEIGHT OF SIGHT LINE ABOVE BORE CENTERLINE (INCHES);
130 REM GUN ELEVATION ANGLE (MINUTES); RATIO AIR DENSITY TO
140 REM STANDARD DENSITY; AIR TEMPERATURE (DEG F); RANGE PRINT
150 REM INTERVAL (YARDS/METERS); RANGE TO TERMINATE TRAJECTORY (YARDS/METERS);
160 REM RANGE WIND SPEED (MPH--POSITIVE IF WIND BLOWS FROM GUN TO
170 REM TARGET); AND CROSSWIND SPEED (MPH--POSITIVE IF WIND BLOWS
180 REM FROM LEFT TO RIGHT ACROSS FIRING LINE).
190 REM THE PROGRAM ALSO PROVIDES AN OPTION TO ADJUST THE TRAJECTORY
200 REM TO PASS THROUGH A SPECIFIED POINT IN SPACE, DENOTED BY
210 REM RMATCH (YARDS/METERS), AND HMATCH (INCHES).  THE GUN ELEVATION ANGLE
220 REM IS ADJUSTED UNTIL THE TRAJECTORY PASSES THROUGH THE POINT
230 REM (RMATCH,HMATCH).  IF NO ADJUSTMENT IS DESIRED, INPUT ZEROS
240 REM FOR RMATCH AND HMATCH, AND THE TRAJECTORY WILL BE RUN WITH
250 REM THE INPUT ELEVATION ANGLE.
260 REM THE PROGRAM SOLVES THE TRAJECTORY BY NUMERICAL INTEGRATION
270 REM USING THE HEUN METHOD, WHICH IS AN ITERATIVELY APPLIED
280 REM SECOND ORDER PREDICTOR-CORRECTOR TECHNIQUE.
290 REM
300 REM THE TRAJECTORY OUTPUT IS RANGE (YARDS/METERS); HEIGHT RELATIVE TO
310 REM LINE OF SIGHT (INCHES); DEFLECTION (INCHES); TOTAL VELOCITY
320 REM (FT/SEC); TIME OF FLIGHT (SECONDS); RANGE COMPONENT OF
330 REM VELOCITY (VX--FT/SEC); VERTICAL COMPONENT OF VELOCITY
340 REM (VY--FT/SEC); AND HORIZONTAL COMPONENT OF VELOCITY (VZ--FT/SEC).
350 REM
360 REM DEFINE PROGRAM CONSTANTS
370 REM
380 G = 32.174#
390 E1 = .00001#
400 U1$ = "###### #######.#     #####.#   #####.#  ###.###  #####.#  #####.# ####.#"
410 REM INPUT DRAG COEFFICIENT TABLE.
420 REM
430 DIM M(50), D(50), E(21), H(21)
440 DIM Q(101), R(101), T(101), V(101)
450 DIM W(101), X(101), Y(101), Z(101)
460 PRINT
470 PRINT "ENTER DRAG FUNCTION TO BE USED:";
480 INPUT K$
490 PRINT
500 PRINT "THE DRAG COEFFICIENT IS ENTERED AS A"
510 PRINT "TABLE OF MACH NUMBER (M) VERSUS"
520 PRINT "DRAG COEFFICIENT (CD).  ENTER THE DATA"
530 PRINT "AS -- M, CD [RETURN] -- BEGINNING"
540 PRINT "WITH THE LOWEST MACH NUMBER."
550 PRINT
560 PRINT "THE LAST (HIGHEST) MACH NUMBER IN THE TABLE"
570 PRINT "MUST BE ENTERED AS A NEGATIVE VALUE,"
580 PRINT "WHICH SIGNALS THE FINAL DRAG INPUT TO THE PROGRAM."
590 PRINT
600 PRINT
610 I = 1
620 PRINT "ENTER M, CD:";
630 INPUT M(I), D(I)
640 IF M(I) < 0 THEN 670
650 I = I + 1
660 GOTO 620
670 J = I
680 M(J) = ABS(M(J))
690 LPDC = 0
700 CLS
710 PRINT
720 PRINT "DRAG FUNCTION: "; K$
730 PRINT
740 PRINT
750 PRINT "MACH NO.", "  CD"
760 REM ECHO INPUT DRAG COEFFCIENT TABLE.
770 FOR I = 1 TO J
780 PRINT M(I), D(I)
790 NEXT I
800 JJ = J
810 PRINT
820 PRINT
830 PRINT "IS THE ABOVE DRAG TABLE CORRECT?"
840 PRINT "ENTER Y FOR YES, N FOR NO:";
850 INPUT K1$
860 IF K1$ = "N" THEN 460
870 PRINT
880 PRINT
890 PRINT "RANGE IN YARDS OR METERS? (ENTER Y FOR YARDS, M FOR METERS):";
900 INPUT YM$
910 PRINT
920 PRINT
930 PRINT "ARMY STANDARD METRO OR ICAO? (ENTER S FOR STANDARD, I FOR ICAO):";
940 INPUT SI$
950 PRINT
960 PRINT
970 IF YM$ = "Y" THEN 1010
980 UC = 1# / .3048#
990 UR$ = "(METERS)"
1000 GOTO 1030
1010 UC = 3#
1020 UR$ = "(YARDS)"
1030 DINT = 1#
1040 REM DINT IS THE DISTANCE INTEGRATION STEP (YARDS/METERS)
1050 D3 = DINT * UC
1060 IF SI$ = "I" THEN 1140
1070 RH1 = -.00003158#
1080 RH2 = 0#
1090 TK1 = -.000006015#
1100 TK2 = 0#
1110 PIR = -.0002048757#
1120 VV1 = 49.19#
1130 GOTO 1220
1140 RH1 = -.00002926#
1150 RH2 = -.0000000001#
1160 TK1 = -.000006858#
1170 TK2 = -.00000000002776#
1180 PIR = -.000208551#
1190 VV1 = 49.0223#
1200 REM NOTE: PIR = -(PI/8)*(RHO0/144)
1210 REM
1220 PRINT "ENTER PROJECTILE IDENTIFICATION:";
1230 INPUT K2$
1240 PRINT
1250 PRINT "ENTER MUZZLE VELOCITY (FEET/SECOND):";
1260 INPUT V0
1270 PRINT
1280 PRINT "ENTER BALLISTIC COEFFICIENT (LB/IN 2):";
1290 INPUT C
1300 PRINT
1310 PRINT "ENTER HEIGHT OF SIGHT LINE ABOVE BORE LINE (INCHES):";
1320 INPUT H0
1330 PRINT
1340 PRINT "ENTER GUN ELEVATION ANGLE (MINUTES):";
1350 INPUT P0
1360 PRINT
1370 PRINT "ENTER RATIO OF AIR DENSITY TO SEA LEVEL STANDARD:";
1380 INPUT R0
1390 PRINT
1400 PRINT "ENTER AIR TEMPERATURE (DEGREES, FAHRENHEIT):";
1410 INPUT TDF
1420 PT0 = INT(10 * TDF) / 10
1430 PRINT
1440 PRINT "ENTER RANGE PRINT INTERVAL"; UR$; ":";
1450 INPUT P2
1460 PRINT
1470 PRINT "ENTER RANGE TO TERMINATE TRAJECTORY "; UR$; ":";
1480 INPUT R3
1490 PRINT
1500 N1 = INT(R3 / P2 + 1.5)
1510 IF N1 <= 101 THEN 1570
1520 PRINT "THIS PRINT INTERVAL GIVES OVER 100 LINES OF OUTPUT!"
1530 PRINT "INCREASE PRINT STEP? (ENTER Y FOR YES, N FOR NO:)";
1540 INPUT PPI$
1550 PRINT
1560 IF PPI$ = "Y" THEN 1440
1570 PRINT "ENTER RANGE WIND SPEED (MILES/HOUR):";
1580 INPUT W1
1590 PRINT
1600 PRINT "ENTER CROSS WIND SPEED (MILES/HOUR):";
1610 INPUT W3
1620 PRINT
1630 PRINT "THE FOLLOWING TWO INPUTS SPECIFY THE COORDINATES"
1640 PRINT "OF A POINT THROUGH WHICH THE TRAJECTORY MUST PASS"
1650 PRINT "IF THIS OPTION IS NOT DESIRED, INPUT ZERO FOR BOTH"
1660 PRINT "THE MATCH RANGE AND THE MATCH HEIGHT."
1670 PRINT
1680 PRINT "ENTER THE TRAJECTORY MATCH RANGE, RMATCH "; UR$; ":";
1690 INPUT R8
1700 PRINT
1710 PRINT "ENTER THE TRAJECTORY MATCH HEIGHT, HMATCH (INCHES):";
1720 INPUT H8
1730 PRINT
1740 PRINT
1750 REM INITIALIZE INPUT VALUES FOR TRAJECTORY CALCULATION.
1760 N = 2
1770 J = 0
1780 P1 = W1
1790 P3 = W3
1800 P4 = H0
1810 H0 = -(H0)
1820 R4 = D3 * P2
1830 R5 = D3 * R3
1840 R6 = D3 * R8
1850 PR4 = P2
1860 PR5 = R3
1870 PR6 = R8
1880 H3 = H8 / 12
1890 W1 = (22 * W1) / 15
1900 W3 = (22 * W3) / 15
1910 V3 = V0 * COS(P0 / 3437.74677#)
1920 V4 = V0 * SIN(P0 / 3437.74677#)
1930 V5 = 0
1940 R1 = 0
1950 PR1 = 0
1960 H1 = H0 / 12
1970 D1 = 0
1980 T1 = 0
1990 PX3 = PR4
2000 REM CHECK FOR TRAJECTORY MATCH VALUES
2010 IF PR6 <> 0 THEN 2050
2020 L = 0
2030 PR7 = PR5
2040 GOTO 2070
2050 L = 1
2060 PR7 = PR6
2070 IF L = 1 THEN 2470
2080 REM PRINT HEADERS AND FIRST LINE OUTPUT
2090 CLS
2100 PRINT
2110 IF SI$ = "S" THEN PRINT "ARMY STANDARD METRO"
2120 IF SI$ = "I" THEN PRINT "ICAO STANDARD ATMOSPHERE"
2130 PRINT
2140 PRINT "DRAG FUNCTION:  "; K$
2150 PRINT
2160 PRINT "PROJECTILE IDENTIFICATION:  "; K2$
2170 PRINT
2180 PRINT "MUZ VEL", "   C", "   H0", "  ELEV", "DENSITY"
2190 PRINT "(FT/SEC)", "(LB/IN2)", "(INCHES)", "(MINUTES)", " RATIO"
2200 PRINT
2210 PP0 = INT(1000 * P0 + .5) / 1000
2220 PRINT V0, C, P4, PP0, R0
2230 PRINT
2240 PRINT
2250 PRINT " TEMP", "RANGEWIND", "CROSSWIND", "RMATCH", "HMATCH"
2260 PRINT "(DEG,F)", "  (MPH)", "  (MPH)", UR$, "(INCHES)"
2270 PRINT
2280 PRINT PT0, P1, P3, R8, H8
2290 PRINT
2300 PRINT
2310 PRINT
2320 PRINT TAB(1); "RANGE"; TAB(11); "HEIGHT"; TAB(23); "DEFL."; TAB(33); "VEL"; TAB(41); "TIME"; TAB(51); "VX"; TAB(60); "VY"; TAB(68); "VZ"
2330 PRINT UR$; TAB(12); "(IN)"; TAB(24); "(IN)"; TAB(32); "(FPS)"; TAB(40); "(SEC)"; TAB(50); "(FPS)"; TAB(59); "(FPS)"; TAB(67); "(FPS)"
2340 PRINT
2350 Q(1) = R1
2360 R(1) = H0
2370 T(1) = D1
2380 V(1) = V0
2390 W(1) = T1
2400 X(1) = V3
2410 Y(1) = V4
2420 Z(1) = V5
2430 PRINT USING U1$; Q(1); R(1); T(1); V(1); W(1); X(1); Y(1); Z(1)
2440 REM
2450 REM BEGIN TRAJECTORY CALCULATION
2460 REM
2470 C3 = (PIR * R0) / C
2480 B1 = SQR((V3 - W1) ^ 2 + V4 ^ 2 + (V5 - W3) ^ 2)
2490 T0 = (TDF + 459.67#) * EXP((TK1 + TK2 * H1) * H1) - 459.67#
2500 V1 = VV1 * SQR(T0 + 459.67#)
2510 X1 = B1 / V1
2520 GOSUB 3990
2530 C1 = X2
2540 C4 = (C3 * C1 * B1 * EXP((RH1 + RH2 * H1) * H1)) / V3
2550 A1 = C4 * (V3 - W1)
2560 A2 = C4 * V4 - G / V3
2570 A3 = C4 * (V5 - W3)
2580 REM APPLY EULER PREDICTOR FORMULA
2590 R2 = R1 + D3
2600 PR2 = PR1 + DINT
2610 V6 = V3 + A1 * D3
2620 V7 = V4 + A2 * D3
2630 V8 = V5 + A3 * D3
2640 B2 = SQR((V6 - W1) ^ 2 + V7 ^ 2 + (V8 - W3) ^ 2)
2650 REM APPLY ITERATIVE HEUN CORRECTOR FORMULA
2660 U = B2
2670 T0 = (TDF + 459.67#) * EXP((TK1 + TK2 * H1) * H1) - 459.67#
2680 V1 = VV1 * SQR(T0 + 459.67#)
2690 X1 = B2 / V1
2700 GOSUB 3990
2710 C2 = X2
2720 C5 = (C3 * C2 * B2 * EXP((RH1 + RH2 * H1) * H1)) / V6
2730 A4 = C5 * (V6 - W1)
2740 A5 = C5 * V7 - G / V6
2750 A6 = C5 * (V8 - W3)
2760 V6 = V3 + .5 * (A1 + A4) * D3
2770 V7 = V4 + .5 * (A2 + A5) * D3
2780 V8 = V5 + .5 * (A3 + A6) * D3
2790 B2 = SQR((V6 - W1) ^ 2 + V7 ^ 2 + (V8 - W3) ^ 2)
2800 E2 = ABS((B2 - U) / B2)
2810 IF E2 > E1 THEN 2660
2820 REM COMPUTE VALUES AT R2.
2830 H2 = H1 + ((V4 + V7) / (V3 + V6)) * D3
2840 D2 = D1 + ((V5 + V8) / (V3 + V6)) * D3
2850 T2 = T1 + (2 * D3) / (V3 + V6)
2860 V2 = SQR(V6 ^ 2 + V7 ^ 2 + V8 ^ 2)
2870 REM RESET CONDITIONS AT R1 TO NEW CONDITIONS AT R2.
2880 R1 = R2
2890 PR1 = PR2
2900 H1 = H2
2910 D1 = D2
2920 T1 = T2
2930 V3 = V6
2940 V4 = V7
2950 V5 = V8
2960 A1 = A4
2970 A2 = A5
2980 A3 = A6
2990 P5 = PR2
3000 P6 = 12 * H2
3010 P7 = 12 * D2
3020 REM CHECK STATUS OF PRINT CONDITIONS.
3030 IF L = 1 THEN 3180
3040 IF PR2 < PX3 THEN 3180
3050 Q(N) = P5
3060 R(N) = P6
3070 T(N) = P7
3080 V(N) = V2
3090 W(N) = T2
3100 X(N) = V6
3110 Y(N) = V7
3120 Z(N) = V8
3130 PRINT USING U1$; Q(N); R(N); T(N); V(N); W(N); X(N); Y(N); Z(N)
3140 N = N + 1
3150 REM INCREMENT PRINT RANGE.
3160 PX3 = PX3 + PR4
3170 REM CHECK FOR CONDITIONS TO STOP TRAJECTORY CALCULATION
3180 IF PR2 >= PR7 THEN 3200
3190 GOTO 2590
3200 IF L = 0 THEN 3400
3210 REM INTERATION TO ADJUST ELEVATION FOR MATCH RANGE AND HEIGHT.
3220 E3 = ABS(H2 - H3)
3230 IF E3 < .00001# THEN 3350
3240 J = J + 1
3250 E(J) = P0
3260 H(J) = H2
3270 IF J >= 2 THEN 3320
3280 E(J + 1) = P0 + 2#
3290 P0 = E(J + 1)
3300 GOTO 1910
3310 REM ALGORITHM TO ADJUST ELEVATION ANGLE
3320 E(J + 1) = E(J) + (H3 - H(J)) * (E(J - 1) - E(J)) / (H(J - 1) - H(J))
3330 IF J < 20 THEN 3290
3340 REM FINAL TRAJECTORY WITH PRINTS
3350 PR6 = 0
3360 GOTO 1910
3370 REM [CHECK FOR HARD COPY OF OUTPUT]
3380 REM [SAVE OUTPUT FOR HARD COPY, IF DESIRED]
3390 REM (N2=NO. OF TRAJECTORY LINES IN OUTPUT)
3400 N2 = N - 1
3410 PRINT
3420 PRINT
3430 PRINT "DO YOU WANT HARD COPY OF THIS OUTPUT?"
3440 PRINT "ENTER Y FOR YES, N FOR NO:";
3450 INPUT K1$
3460 IF K1$ = "N" THEN 3880
3470 IF K1$ = "Y" THEN 3490
3480 GOTO 3420
3490 REM ECHO INPUT DRAG COEFFICIENT TABLE.
3500 LPRINT
3510 LPRINT "DRAG FUNCTION:  "; K$
3520 LPRINT
3530 IF LPDC = 1 THEN 3590
3540 LPRINT "MACH NO.", "  CD"
3550 FOR I = 1 TO JJ
3560 LPRINT M(I), D(I)
3570 NEXT I
3580 LPDC = 1
3590 LPRINT
3600 IF SI$ = "S" THEN LPRINT "ARMY STANDARD METRO"
3610 IF SI$ = "I" THEN LPRINT "ICAO STANDARD ATMOSPHERE"
3620 LPRINT
3630 LPRINT "PROJECTILE IDENTIFICATION:  "; K2$
3640 LPRINT
3650 LPRINT "MUZ VEL", "   C", "   H0", "  ELEV", "DENSITY"
3660 LPRINT "(FT/SEC)", "(LB/IN2)", "(INCHES)", "(MINUTES)", " RATIO"
3670 LPRINT
3680 PP0 = INT(1000 * P0 + .5) / 1000
3690 LPRINT V0, C, P4, PP0, R0
3700 LPRINT
3710 LPRINT
3720 LPRINT " TEMP", "RANGEWIND", "CROSSWIND", "RMATCH", "HMATCH"
3730 LPRINT "(DEG,F)", "  (MPH)", "  (MPH)", UR$, "(INCHES)"
3740 LPRINT
3750 LPRINT PT0, P1, P3, R8, H8
3760 LPRINT
3770 LPRINT
3780 LPRINT TAB(1); "RANGE"; TAB(11); "HEIGHT"; TAB(23); "DEFL."; TAB(33); "VEL"; TAB(41); "TIME"; TAB(51); "VX"; TAB(60); "VY"; TAB(68); "VZ"
3790 LPRINT UR$; TAB(12); "(IN)"; TAB(24); "(IN)"; TAB(32); "(FPS)"; TAB(40); "(SEC)"; TAB(50); "(FPS)"; TAB(59); "(FPS)"; TAB(67); "(FPS)"
3800 LPRINT
3810 FOR N = 1 TO N2
3820 LPRINT USING U1$; Q(N); R(N); T(N); V(N); W(N); X(N); Y(N); Z(N)
3830 NEXT N
3840 REM CHECK FOR ADDITIONAL CASES.
3850 LPRINT
3860 LPRINT
3870 LPRINT
3880 PRINT
3890 PRINT
3900 PRINT "RUN ANOTHER CASE? (ENTER Y FOR YES, N FOR NO:)";
3910 INPUT K1$
3920 IF K1$ = "N" THEN 4070
3930 PRINT
3940 PRINT "USE SAME DRAG COEFFICIENT TABLE?"
3950 PRINT "ENTER Y FOR YES, N FOR NO:";
3960 INPUT K1$
3970 IF K1$ = "N" THEN 460
3980 GOTO 870
3990 REM SUBROUTINE FOR INTERPOLATION IN DRAG TABLE
4000 I = 1
4010 IF X1 < M(I + 1) THEN 4040
4020 I = I + 1
4030 GOTO 4010
4040 S = (D(I + 1) - D(I)) / (M(I + 1) - M(I))
4050 X2 = D(I) + S * (X1 - M(I))
4060 RETURN
4070 END