100 REM GALAXY
101 REM
102 REM This program is listed in Appendix J of "An Introduction to
103 REM Modern Astrophysics," Bradley W. Carroll and Dale A. Ostlie,
104 REM Addison-Wesley Publishing Company, copyright 1996.
105 REM
110 REM A galaxy collision program adapted from the listing
120 REM in Astronomy magazine, Dec 1988, by M. C. Schroeder
130 REM Neil F. Comins.  This program allows two galaxies
140 REM to interact via their gravitational fields.  The 
150 REM nuclei of the TARGET galaxy and the INTRUDER galaxy
160 REM are point masses, and they move according to their
170 REM gravitational attraction.  The TARGET galaxy has a
180 REM disk of stars.  The stars respond to the gravity of
190 REM the nuclei but do not influence their motions.  The
200 REM dimension statements below allow up to 1000 stars.
210 REM You can use more by redimensioning these arrays.
220 REM
230 DIM X(1000), Y(1000), VX(1000), VY(1000), Z(1000), VZ(1000)
240 CLS
250 REM
260 REM The stars are distributed
270 REM in the TARGET galaxy (only) in rings
280 REM The total number of stars
290 REM is the product of the \# of rings and
300 REM the \# of stars per ring.
310 REM
320 LOCATE 1: PRINT "Input the \# of rings of stars ";
330 PRINT "in the TARGET galaxy (at least 2):  ";
340 INPUT "", NRR
350 LOCATE 5: PRINT "Input the \# of stars per ring ";
360 PRINT "in the TARGET galaxy:  ";
370 INPUT "", NRS
380 NS = NRR * NRS: DR = 20 / (NRR - 1)
390 LOCATE 9: PRINT "Input mass fraction ";
400 PRINT "of INTRUDER galaxy in units of ";
410 PRINT "the target galaxy mass: ";
420 INPUT "", M2
430 LOCATE 13: PRINT "Input the initial X,Y,Z ";
440 PRINT "coordinates of the INTRUDER galaxy:"
450 PRINT "The TARGET galaxy is located at 0,0,0.  ";
460 INPUT "", X2, Y2, Z2
470 LOCATE 17: PRINT "Input the initial X,Y,Z ";
480 PRINT "velocities of the INTRUDER galaxy:"
490 PRINT "The TARGET galaxy is initially at rest.  ";
500 INPUT "", VX2, VY2, VZ2
510 REM The position of the TARGET galaxy
520 REM is assigned by the computer.
530 REM
540 LOCATE 21: PRINT "Input \# of time steps ";
550 PRINT "for this run:  ";
560 INPUT "", NTSPR
570 LOCATE 24: PRINT "Are these inputs correct? ";
580 PRINT "(1=yes,0=no):  ";
590 INPUT "", ANS
600 IF ANS <> 1 THEN 240
610 SCREEN 1, 0, 0
620 WINDOW (0, 0)-(300, 200)
630 LOCATE 12: PRINT "        Thinking...."
640 REM Initialize TARGET galaxy mass,
650 REM position, and velocity.
660 M1 = 5: X1 = 150: Y1 = 100: VX1 = 0: VY1 = 0: Z1 = 0: VZ1 = 0: SF2 = 2
670 REM Scale the INTRUDER galaxy mass
680 REM position and velocities.
690 TIME = 0: M2 = M2 * M1: X2 = X2 + X1: Y2 = Y2 + Y1: Z2 = Z2 + Z1
700 REM Set up initial load,
710 REM place stars in concentric rings
720 REM from r=10 to r=30 at
730 REM intervals of dr.
740 FOR IR = 1 TO NRR
750 R = 10 + (IR - 1) * DR
760 V = SQR(M1 / R): TH = (.5 * V / R) * (180 / 3.14159)
770 IF R = 10 THEN V = .9 * V
780 FOR IT = 1 TO NRS
790 T = (IT - 1) * 360 / NRS
800 T1 = 3.14159 * (T - TH) / 180
810 I = I + 1
820 REM Initialize star positions.
830 X(I) = R * COS(T / 57.2958) + 150
840 Y(I) = R * SIN(T / 57.2958) + 100
850 VZ(I) = 0: Z(I) = 0
860 REM Initialize star velocities
870 REM to place them in circular
880 REM orbits about the target galaxy.
890 VX(I) = -V * SIN(T1)
900 VY(I) = V * COS(T1)
910 NEXT IT
920 NEXT IR
930 GOSUB 1440
940 REM Particle pusher routine.
950 FOR K = 1 TO NTSPR
960 FOR J = 1 TO 1
970 FOR I = 1 TO NS
980 REM
990 REM Determine the force on a star
1000 REM due to the galactic centers.
1010 REM SF2 below, called the softening factor,
1020 REM is used to prevent overflows
1030 REM during force calculations.
1040 REM SF2 is assigned a value above.
1050 REM You may experiment with its value
1060 REM but it should be kept as small as possible
1070 REM to better approximate
1080 REM a true 1/r**2 force.
1090 REM
1100 R1 = M1 / ((X(I) - X1) ^ 2 + (Y(I) - Y1) ^ 2 + (Z(I) - Z1) ^ 2 + SF2) ^ 1.5
1110 R2 = M2 / ((X(I) - X2) ^ 2 + (Y(I) - Y2) ^ 2 + (Z(I) - Z2) ^ 2 + SF2) ^ 1.5
1120 REM Calculate star's x,y,z accelerations.
1130 AX = R1 * (X1 - X(I)) + R2 * (X2 - X(I))
1140 AY = R1 * (Y1 - Y(I)) + R2 * (Y2 - Y(I))
1150 AZ = R1 * (Z1 - Z(I)) + R2 * (Z2 - Z(I))
1160 REM Update star positions and velocities
1170 REM using a time-centered
1180 REM leap-frog algorithm.
1190 VX(I) = VX(I) + AX
1200 VY(I) = VY(I) + AY
1210 VZ(I) = VZ(I) + AZ
1220 X(I) = X(I) + VX(I)
1230 Y(I) = Y(I) + VY(I)
1240 Z(I) = Z(I) + VZ(I)
1250 NEXT I
1260 REM Update positions
1270 REM and velocities of the galactic centers.
1280 S = ((X1 - X2) ^ 2 + (Y1 - Y2) ^ 2 + (Z1 - Z2) ^ 2 + SF2) ^ 1.5
1290 AX = (X2 - X1) / S: AY = (Y2 - Y1) / S: AZ = (Z2 - Z1) / S
1300 VX1 = VX1 + M2 * AX: VY1 = VY1 + M2 * AY: VZ1 = VZ1 + M2 * AZ
1310 VX2 = VX2 - M1 * AX: VY2 = VY2 - M1 * AY: VZ2 = VZ2 - M1 * AZ
1320 X1 = X1 + VX1
1330 Y1 = Y1 + VY1: Z1 = Z1 + VZ1
1340 X2 = X2 + VX2
1350 Y2 = Y2 + VY2: Z2 = Z2 + VZ2
1360 TIME = TIME + 1
1370 NEXT J
1380 GOSUB 1440
1390 NEXT K
1400 LOCATE 23: PRINT "Continue (1=yes,0=no)"; : INPUT "", ANS
1410 IF ANS = 1 THEN 950
1420 GOTO 1730
1430 REM Update screen display.
1440 CLS
1450 REM Calculate center of mass of system
1460 REM for use as center
1470 REM of output display.
1480 XC = (M1 * X1 + M2 * X2) / (M1 + M2)
1490 YC = (M1 * Y1 + M2 * Y2) / (M1 + M2)
1500 ZC = (M1 * Z1 + M2 * Z2) / (M1 + M2)
1510 REM Set up output display.
1520 LINE (0, 50)-(300, 50)
1530 LINE (150, 200)-(150, 50)
1540 REM Calculate position of galactic centers
1550 REM and display on screen.
1560 XX1 = (X1 - XC): YY1 = (Y1 - YC): ZZ1 = (Z1 - ZC)
1570 XX2 = (X2 - XC): YY2 = (Y2 - YC): ZZ2 = (Z2 - ZC)
1580 CIRCLE ((XX1 + 75), (YY1 + 125)), 2, 1
1590 CIRCLE ((XX1 + 225), (2 * ZZ1 + 125)), 2, 2
1600 CIRCLE ((XX2 + 75), (YY2 + 125)), 1, 3
1610 CIRCLE ((XX2 + 225), (2 * ZZ2 + 125)), 1, 3
1620 REM Place stars on screen.
1630 FOR I = 1 TO NS
1640 XP = (X(I) - XC): YP = (Y(I) - YC): ZP = 2 * (Z(I) - ZC)
1650 PSET ((XP + 75), (YP + 125)), 2: PSET ((XP + 225), (ZP + 125)), 1
1660 NEXT I
1670 REM Print out display headings.
1680 LOCATE 21: PRINT TAB(2);
1690 PRINT "       x-y", "               x-z"
1700 RT = 1.2 * TIME
1710 LOCATE 23: PRINT "step="; TIME; "  time="; RT; "million years";
1720 RETURN
1730 WIDTH 80: SCREEN 0: CLS : END