10 ' TRIBODY.BAS R.S. FRITZIUS 05/05/1990 revised 24-JAN-1999 13 Jun 2009 'IDELAY, IDELMIN, IDELMAX, AND J were changed to 'DELAY!, DELMIN!, DELMAX!, AND DJ! on 2 June 2007 cls:print:print print " ΙΝΝΝΝΝΝΝΝΝ» " print " ΘΝΝΝ» ΙΝΝΝΌ " print " ΙΝ» ΙΝΝΝΝΝΝΝΝΝ» ΙΝ» ΙΝ» Ί Ί ΙΝΝΝΝΝΝΝΝ» ΙΝΝΝΝΝΝΝΝ» " print " Ί Ί Ί ΙΝΝΝΝΝ» Ί Ί Ί Ί Ί Ί Ί Ί ΙΝΝΝΝ» Ί Ί ΙΝΝΝΝΝΝΌ " print " Ί Ί Ί Ί Ί Ί Ί Ί ΙΝ» Ί Ί ΙΝΝ» Ί Ί Ί ΘΝΝΝΝΌ Ί Ί Ί " print " Ί Ί Ί Ί Ί Ί Ί Ί Ί Ί Ί Ί ΘΝΝΌ Ί Ί Ί ΙΝΝΝΝΝΝΌ Ί Ί " print " Ί Ί Ί ΘΝΝΝΝΝΌ Ί Ί ΘΝΌ ΘΝΌ Ί Ί Ί Ί ΘΝΝΝΝΝΝ» Ί ΘΝΝΝΝΝΝ» " print " Ί Ί ΘΝΝΝΝΝΝΝΝΝΌ ΘΝΝΝΝΝΝΝΝΝΌ ΘΝΌ ΘΝΝΝΝΝΝΝΝΌ ΘΝΝΝΝΝΝΝΝΌ " print " Ί ΘΝΝΝΝΝΝ» TM " print " ΘΝΝΝΝΝΝΝΝΌ Scientific Simulations Software" print print " T R I B O D Y" 'on error goto 10000 DEFINT I,J,K,N DX=0 ' Delta X between pairs of bodies DY=0 ' Delta Y between pairs of bodies I=0 ' General Counter J=0 ' General Counter K=0 ' General Counter C1=20 C2=-9 DT=1 ' Time Increment for position updates FI=0 ' Force on Ith element ICOUNT=0 ' Counter used for dashed line center of mass plot DIM X(3),Y(3) ' X and Y coordinate arrays for bodies DIM M(3) ' Mass of bodies array DIM R(3,3) ' temporary DIM VX(3),VY(3) ' X and Y velocity arrays for bodies DELAY!=10 ' Dummy value DELMIN!=2 'Delay Minimum DELMAX!=25000'Delay Maximum IL=21 ' Input Line GOSUB 3200 ' Calculate operational delay ESC$=chr$(27) ' Escape Key 20 CLS FOR I=1 TO 200 ' Pause before printing next screen FOR DJ!=1 TO DELAY!:NEXT DJ! NEXT I PRINT:PRINT:PRINT PRINT " NEWTONIAN THREE-BODY MODELING" PRINT PRINT " (TRIBODY)" PRINT 30 PRINT " Copyright (c) 1992 Robert S. Fritzius" PRINT 40 PRINT " Portions of run-time module Copyrighted by Microsoft, 1982-1987" PRINT PRINT " Manufactured by" PRINT PRINT " Low-Tec 305 Hillside Drive Starkville, MS 39759" PRINT print print print PRINT 70 locate 23,18:print "S=Simulations I=Info ESC=Exit program" 80 Z$=INKEY$:IF Z$="" THEN 80 85 if z$="S" or z$="s" then 100 87 if z$="I" or z$="i" then gosub 5000:goto 100 90 if z$=ESC$ then cls:SCREEN 0:END 92 goto 80 100 ' Statement 100 is the initial conditions starting point FI=0 ' G*Mj / R^2 actually (G*Mi*Mj)/(Mi*R^2) = dV/dt I1=0 ' Counter in operational delay calculation routine I2=0 ' Counter in operational delay calculation routine Ai=0 ' Acceleration of Ith body KC12=0 ' Center of Mass (bodies 1 & 2) plot ON/OFF 1=ON 0=OFF M12=0 ' Mass of bodies 1 and 2 combined MALL=0 ' Mass of all bodies combined NB=3 ' Number of Bodies NP=0 ' Number of Bodies to plot with respect to C.M. 1,2 230 R=0 ' Distance between pairs of bodies 240 RR=0 ' Distance Squared - Body to Body in arbitrary units R1=0 ' Radius body 1 R2=0 ' Radius body 2 R3=0 ' Radius body 3 250 S1=1 ' Scale Factor 1 (General X,Y) 255 S2=1 ' Scale Factor 2 (for center of mass plot) 260 SF=0 ' This will be the adjustable Scale Factor 265 SN=0 ' Senario Number 515 TB=20 ' Tab Setting on Menu Screen 270 TITLE1$="" ' Scenario Title-1" 275 TITLE2$="" ' Scenario Title-2 277 COM12$="" ' Center of Gravity bodies 1 & 2 label 280 IXREF=310 ' X Origin X/Y Plot Screen Laboratory coordinates 290 JYREF=95 ' Y Origin X/Y Plot screen Laboratory coordinates 300 IXQ=100 ' X Origin X/Y plot for motion about c.m. 1,2 310 JYQ=100 ' Y Origin X/Y plot for motion about c.m. 1,2 320 IXS=0 ' X on Screen (Integer) 330 JYS=0 ' Y on Screen (Integer) 340 X12=0 ' X center of mass bodies 1,2 350 Y12=0 ' Y center of mass bodies 1,2 360 XALL=0 ' X center of mass all bodies 370 YALL=0 ' Y center of mass all bodies 375 ILB=20 ' X coordinate for C.O.G. Label for bodies 1 & 2 RMIN=0 ' Minimum radius of planet or star ILIMIT=30000 410 UA$=CHR$(0)+CHR$(72) ' Up Arrow 420 DA$=CHR$(0)+CHR$(80) ' Down Arrow 430 LA$=CHR$(0)+CHR$(75) ' Left Arrow 440 RA$=CHR$(0)+CHR$(77) ' Right Arrow 445 F1$=CHR$(0)+CHR$(59) ' F1 Key 450 F2$=CHR$(0)+CHR$(60) ' F2 Key 455 F3$=CHR$(0)+CHR$(61) ' F3 Key 460 F4$=CHR$(0)+CHR$(62) ' F4 Key 465 F9$=CHR$(0)+CHR$(67) ' F9 Key 470 SP$=CHR$(32) ' Space Bar 480 IXMAX=680 490 JYMAX=210 500 DIM VN(15) ' VN() array is used for changing variable parameters 'DIM IXS(3),JYS(3) 502 BL$=" Space Bar to Stop/Start/Change " 504 BLANK$=STRING$(78," ") 510 CLS:SCREEN 0 '-----------menu screen----------- 520 LOCATE 5,1 PRINT TAB(31) "M E N U":PRINT PRINT TAB(TB) "a. Binary Star and Interloper" PRINT TAB(TB) "b. Sun, Planet, and Moon" PRINT TAB(TB) "c. Sun, Planet, and Comet" PRINT TAB(TB) "d. Info PRINT TAB(TB) " ESC Exit program" LOCATE 13,TB:PRINT " Which selection ? "; 580 WS$=INKEY$:IF WS$="" THEN 580 IF WS$=ESC$ THEN CLS:SCREEN 0:END WS%=ASC(WS$) IF WS%<64 OR WS%>122 THEN 580 'Outside usable range IF WS%>96 AND WS%<123 THEN WS%=WS%-32 ' make upper case WS%=WS%-64 ' Downshift 590 IF WS%<1 OR WS%>5 THEN 580 597 SN=0 'Selection no.s are set in scenario setup sections 600 ON WS% GOTO 620,740,860,605 605 gosub 5000:goto 510 620 '-------scenario 1------------ SN=1 TITLE1$="INTERLOPER" TITLE2$=" STAR-1 STAR-2 INTERLOPER" ' TITLE2$=" PLANET MOON SUN" ' TITLE2$=" SUN PLANET COMET" ' "Initial X (a) ###.### (b) ###.### (c) ###.### " COM12$="Star-1 & Star-2" S1=2:S2=3 NP=2 670 KC12=1 ' Center of mass plot ON 672 IXEL=-2:IXER=7 674 JYEB=-4:JYET=4 690 STP=1 X1= 0 :X2= 0 :X3= 6 Y1= 1.5 :Y2= -2.3 :Y3= 72 M1= 2 :M2= 1 :M3= 0.2 U1= 0.03 :U2= -0.05 :U3= 0 V1= 0 :V2= 0 :V3= -0.2 G=.02 RMIN=.001 722 IXREF=320:JYREF=100 'regular plot lower limits 725 IXQ=100:JYQ=75 'c.o.m. 730 GOTO 980 740 '-------scenario 2 ------------ 745 SN=2 750 TITLE1$="PLANET, MOON, and SUN" 760 TITLE2$=" PLANET MOON SUN" 765 COM12$="Planet & Moon" 770 S1=.5:S2=4 785 NP=2 790 KC12=0 ' Center of mass plot OFF 792 IXEL=-10:IXER=18 794 JYEB=-10:JYET=14 796 STP=4 'planet moon sun X1= 3 :X2= 3 :X3= 3 Y1=-9.1 :Y2=-10.1 :Y3= 3.5 M1= 1 :M2= 0.1 :M3= 25 U1= 0.125:U2= 0.222:U3= -0.005 V1= 0 :V2= 0 :V3= 0 G=.01 RMIN=.001 842 IXREF=380:JYREF=100 'lower limits main xy plot 845 IXQ=100:JYQ=90 'c.o.m. 850 GOTO 980 860 '--------scenario 3------------- 865 SN=3 870 TITLE1$="SUN, PLANET, and COMET" 880 TITLE2$=" SUN PLANET COMET" 885 COM12$="Sun, Planet & Comet" 890 S1=2:S2=5 905 NP=3 910 KC12=1 ' Center of mass plot ON 912 IXEL=-3:IXER=4 914 JYEB=-4:JYET=4 916 STP=1 RMIN=.001 918 'sun earth comet X1= 0 :X2= 0 :X3= 0 Y1= -2 :Y2= -1 :Y3= 3.7 M1= 400 :M2= .1 :M3= .02 R1=.01 :R2= .0001:R3= .00001 U1= 0 :U2=-.1 :U3= -.020 V1= 0 :V2= 0 :V3= 0 G=.000025 962 IXREF=400:JYREF=100 ' lower limits main x,y plot 965 IXQ=130:JYQ=100 ' c.o.m. 970 : 980 REM Load initial variables into parameter changing arrays 985 VN( 1)=X1: VN( 2)=X2: VN( 3)=X3 ' X components 990 VN( 4)=Y1: VN( 5)=Y2: VN( 6)=Y3 ' Y components 1000 VN( 7)=M1: VN( 8)=M2: VN( 9)=M3 ' Masses 1010 VN(10)=U1: VN(11)=U2: VN(12)=U3 ' X components of velocities 1020 VN(13)=V1: VN(14)=V2: VN(15)=V3 ' Y components of velocities ' 1040 CLS:SCREEN 0 SF=S1 LOCATE 3,1 1100 PRINT "Initial parameters for "TITLE1$":" 1110 PRINT 1130 PRINT TITLE2$ 1140 PRINT 1150 PRINT USING "Initial X (a) ###.### (b) ###.### (c) ###.### ";VN( 1);VN( 2);VN( 3) PRINT 1170 PRINT USING "Initial Y (d) ###.### (e) ###.### (f) ###.### ";VN( 4);VN( 5);VN( 6) PRINT 1190 PRINT USING "Mass (g) ###.### (h) ###.### (i) ###.### ";VN( 7);VN( 8);VN( 9) PRINT 1210 PRINT USING "Initial VX (j) ###.### (k) ###.### (l) ###.### ";VN(10);VN(11);VN(12) PRINT 1230 PRINT USING "Initial VY (m) ###.### (n) ###.### (o) ###.### ";VN(13);VN(14);VN(15) PRINT:PRINT 1252 PRINT "ESC Return to main menu (q) Info (r) Restore Demo Settings "; 1280 LOCATE IL,1 PRINT "Which option or variable to change "; 1285 WI$=INKEY$:IF WI$="" THEN 1285 IF WI$=ESC$ THEN 510 WI%=ASC(WI$) IF WI%=13 THEN GOTO 1430 ' ENTER was depressed - run simulation IF WI%<64 OR WI%>122 THEN 1285 ' Outside range of acceptible keystrokes IF WI%>96 AND WI%<123 THEN WI%=WI%-32 'Make it upper case WI%=WI%-64 ' Downshift 1300 IF WI%>18 THEN GOTO 100 IF WI%=16 THEN 1040 ' Don't let Grav Const be a selectable item 1345 IF WI%=17 THEN GOSUB 5900:GOTO 1040 ' (q) Info then to the top at return IF WI%=18 THEN ON SN GOTO 620,740,860 '(r) Restore initial conditions LCWI$=CHR$(WI%+96) 'for screen display of selected item 1355 LOCATE IL,1:PRINT BLANK$; LOCATE IL,1:PRINT " New value for ("LCWI$") "; 1375 INPUT VN$:IF VN$="" THEN 1280 1380 VN=VAL(VN$) LOCATE IL,1:PRINT BLANK$; IF WI%>0 AND WI%<16 THEN VN(WI%)=VN: GOTO 1040 1390 'IF WI%=16 AND VN>0 AND VN<=1 THEN G=VN:GOTO 1040 1410 GOTO 1040 1420 : 1430 CLS:SCREEN 2 1440 REM Load variables into operational arrays 1445 FOR I=1 TO 3 1450 X(I)=VN(I+ 0) 1460 Y(I)=VN(I+ 3) 1470 M(I)=VN(I+ 6) 1480 VX(I)=VN(I+ 9) 1490 VY(I)=VN(I+12) 1495 NEXT I 1500 : 1510 M12=M(1)+M(2) ' Combined Mass of bodies 1 & 2 1550 MALL=M(1)+M(2)+M(3) ' Combined Mass of all three bodies 1560 : 1570 CLS:PRINT " NEWTONIAN THREE BODY INTERACTIONS" ' DRAW "BM"+STR$(IXREF) +","+STR$(JYREF) ' DRAW "M"+STR$(IXREF) +","+STR$(JYREF) ' CIRCLE(IXREF,JYREF),10 1575 'GOSUB 4000 'this plots a grid but takes too much time 1582 ILB=15-LEN(COM12$)/2 1585 LOCATE 21,ILB:PRINT COM12$; 1587 LOCATE 22,4:PRINT "Plotted with respect to"; LOCATE 23,4:PRINT " their Center of Mass"; 1600 LOCATE 24,1 1610 PRINT "ESC=New Run F1=Slower F2=Faster F3=Shrink F4=Expand F9=Erase Arrows=Move Origin"; ' 1670 REM Calculate new X and Y velocities - Loop start point 1675 ICOUNT=ICOUNT+1:IF ICOUNT=40 THEN ICOUNT=0 'see line 1935 FOR I=1 TO NB 'calculate new velocites for Ith mass FOR J=1 TO NB 'based on summation of grav forces of J masses IF I=J THEN 1780 ' No self actions DX=X(I)-X(J):DY=Y(I)-Y(J) RR=DX^2+DY^2 IF RRILIMIT THEN 2700 IF ABS(JYS)>ILIMIT THEN 2700 'IXS(I)=IXS 'JYS(I)=JYS 1860 NEXT I 1870 ' 1880 REM Calculate and plot Center of mass for Bodies 1 and 2 (c.m. 1,2) 1890 X12=(M(1)*X(1)+M(2)*X(2))/M12 1900 Y12=(M(1)*Y(1)+M(2)*Y(2))/M12 1910 IF KC12=1 THEN GOSUB 2300 ' Plot center of mass for bodies 1 & 2 1920 : 1930 REM Calculate and plot Center of mass for all bodies 1935 IF ICOUNT>20 THEN 1980 ' Creates dashed line for C.O.M. all XALL=(M(1)*X(1)+M(2)*X(2)+M(3)*X(3))/MALL YALL=(M(1)*Y(1)+M(2)*Y(2)+M(3)*Y(3))/MALL IF KC12=1 THEN GOSUB 2350 ' Plot center of mass for all bodies 1970 ' 1980 REM Calc and plot positions of bodies 1 & 2 with respect to their c.m. 1990 FOR I=1 TO 2 2000 XC=X(I)-X12 2010 YC=Y(I)-Y12 2030 GOSUB 2400 ' Plot positions 2040 NEXT I GOSUB 2600 ' Print radii 2044 'Crank in some processing delay to keep program 'from going too fast FOR JJ=1 TO 30 FOR DJ!=1 TO DELAY!:NEXT DJ! NEXT JJ 2050 : 2060 LOCATE 2,24:PRINT USING "X(3)=#####.## Y(3)=#####.##";X(3);Y(3); 2062 Z$=INKEY$:IF Z$="" THEN GOTO 1670 IF Z$=ESC$ THEN GOTO 1040 'New Run IF Z$<>SP$ THEN GOTO 2080 2077 Z$=INKEY$:IF Z$="" THEN 2077 'Pause for space bar IF Z$=ESC$ THEN 1040 'New Run 2078 GOTO 1670 2080 IF Z$=UA$ THEN JYREF=JYREF-30: GOTO 1670 2090 IF Z$=DA$ THEN JYREF=JYREF+30: GOTO 1670 2100 IF Z$=LA$ THEN IXREF=IXREF-100:GOTO 1670 2110 IF Z$=RA$ THEN IXREF=IXREF+100:GOTO 1670 2140 IF Z$=F1$ AND DELAY!<0.80*DELMAX! THEN GOSUB 3500:GOTO 1670 ' Slower 2150 IF Z$=F2$ AND DELAY!>1.25*DELMIN! THEN GOSUB 3600:GOTO 1670 ' Faster 2120 IF Z$=F3$ THEN SF=SF/2: GOTO 1670 ' Shrink 2130 IF Z$=F4$ THEN SF=SF*2: GOTO 1670 ' Expand 2155 IF Z$=F9$ THEN GOTO 1570 ' Erase & continue ' 2180 GOTO 1670 ' 2185 '----------------- end of loop --------------------------- ' 2200 REM plot X,Y on screen in graphics format IXS=INT(C1*SF*X(I)+IXREF) JYS=INT(C2*SF*Y(I)+JYREF) IF IXS<0 OR JYS<0 THEN 2260 IF IXS>IXMAX OR JYS>JYMAX THEN 2260 DRAW "BM"+STR$(IXS) +","+STR$(JYS) DRAW " M"+STR$(IXS) +","+STR$(JYS) 2260 RETURN ' 2300 REM plot Center of mass for bodies 1 and 2 IXS=INT(C1*SF*X12+IXREF) JYS=INT(C2*SF*Y12+JYREF) IF IXS<0 OR JYS<0 THEN 2340 IF IXS>IXMAX OR JYS>JYMAX THEN 2340 DRAW "BM"+STR$(IXS) +","+STR$(JYS) DRAW "M"+STR$(IXS) +","+STR$(JYS) 2340 RETURN ' 2350 REM plot Center of mass for all bodies IXS=INT(C1*SF*XALL+IXREF) JYS=INT(C2*SF*YALL+JYREF) IF IXS<0 OR JYS<0 THEN 2390 IF IXS>IXMAX OR JYS>JYMAX THEN 2390 DRAW "BM"+STR$(IXS) +","+STR$(JYS) DRAW "M"+STR$(IXS) +","+STR$(JYS) 2390 RETURN ' 2400 REM plot bodies 1 and 2 with respect to their center of mass IXS=INT(C1*S2*XC + IXQ) JYS=INT(C2*S2*YC + JYQ) IF IXS<0 OR JYS<0 THEN 2470 IF IXS>IXMAX OR JYS>JYMAX THEN 2470 DRAW "BM"+STR$(IXS) +","+STR$(JYS) DRAW "M"+STR$(IXS) +","+STR$(JYS) 2470 RETURN ' 2600 REM Print values of radii IF SN<>3 THEN 2620 IF R(2,3)<.0001 THEN 2620 FREQ=100+1000/R(2,3) SOUND FREQ,1 2620 'LOCATE 3,55:PRINT USING "###.### ";R(1,1);R(1,2);R(1,3); 'LOCATE 4,55:PRINT USING "###.### ";R(2,1);R(2,2);R(2,3); 'LOCATE 5,55:PRINT USING "###.### ";R(3,1);R(3,2);R(3,3); 'LOCATE 6,55:PRINT USING "###,### ";IXS(1);IXS(2);IXS(3); 'LOCATE 7,55:PRINT USING "###,### ";JYS(1);JYS(2);JYS(3); 2660 RETURN ' 2700 'Print overflow message LOCATE 2,10 PRINT "Math overflow is approaching. Press ESC to start new run."; 2710 Z$=INKEY$:IF Z$="" THEN 2710 GOTO 1040 ' 3200 'Calculating processing delay I1=VAL(RIGHT$(TIME$,2)) DJ!=0:K=0 X(1)=10:X(2)=5:Y(1)=20:Y(2)=15 3230 FOR I=1 TO 400 ' dummy math processing DX=X(1)-X(2):DY=Y(1)-Y(2) R=SQR(DX^2+DY^2) NEXT I DJ!=DJ!+1 I2=VAL(RIGHT$(TIME$,2)) IF I2DELMAX! THEN DELAY!=DELMAX! 3310 RETURN 3320 ' 3400 'Print delayed message on screen COLOR 25 ' light blue, flashing locate 17,26:PRINT "Calculating operational speed" COLOR 7 K=1 RETURN 3500 DELAY!=INT(1.25*DELAY!) 'slow down IF DELAY!>DELMAX! THEN DELAY!=DELMAX! 3520 SOUND 1000-INT(40*LOG(DELAY!)),1 3530 RETURN 3590 ' 3600 DELAY!=INT(0.8*DELAY!) ' Speed up IF DELAY!5 then PRINT "Error code ="err:goto 10090 10010 PRINT " This machine does not know how to do high-resolution graphics." PRINT PRINT " Options" PRINT PRINT "(1) Use a PC that has a CGA, VGA, or EGA card. PRINT PRINT "(2) If your PC has a Monochrome Graphics Adapter (MGA) card, run SIMCGA" PRINT " then run this program. (If you don't know what kind of card your" PRINT " machine has, try this option.) When you finish, run RMVCGA, or reset" PRINT " your PC, to remove SIMCGA." PRINT PRINT "SIMCGA, Copyright (C) 1988-1992 by Sydex, Inc., is a program which enables" PRINT "a Monochrome Graphics Adapter (MGA) to simulate a Color Graphics Adapter " PRINT "(CGA). SIMCGA will not work with a Monochrome text Display Adapter (MDA)." PRINT "SIMCGA is licensed to Low-Tec, by Sydex, as a non-supported item. Sydex " PRINT "and Low-Tec may not be held liable for damages arising from the mis-use of" PRINT "SIMCGA." 10090 LOCATE 23,20 PRINT " Press any key to exit program."; 10095 Z$=INKEY$:IF Z$="" THEN 10095 END